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Abstract 


We present observations and timing analyses of 68 millisecond pulsars (MSPs) comprising the 15 yr data set of 
the North American Nanohertz Observatory for Gravitational Waves (NANOGrav). NANOGrav is a pulsar 
timing array (PTA) experiment that is sensitive to low-frequency gravitational waves (GWs). This is 
NANOGrav’s fifth public data release, including both “narrowband” and “wideband” time-of-arrival (TOA) 
measurements and corresponding pulsar timing models. We have added 21 MSPs and extended our timing 
baselines by 3 yr, now spanning nearly 16 yr for some of our sources. The data were collected using the 
Arecibo Observatory, the Green Bank Telescope, and the Very Large Array between frequencies of 327 MHz 
and 3 GHz, with most sources observed approximately monthly. A number of notable methodological and 
procedural changes were made compared to our previous data sets. These improve the overall quality of the 
TOA data set and are part of the transition to new pulsar timing and PTA analysis software packages. For the 
first time, our data products are accompanied by a full suite of software to reproduce data reduction, analysis, 
and results. Our timing models include a variety of newly detected astrometric and binary pulsar parameters, 
including several significant improvements to pulsar mass constraints. We find that the time series of 23 pulsars 
contain detectable levels of red noise, 10 of which are new measurements. In this data set, we find evidence for 
a stochastic GW background. 


Unified Astronomy Thesaurus concepts: Millisecond pulsars (1062); Pulsar timing method (1305); Time series 
analysis (1916); Pulsars (1306); Gravitational waves (678) 


1. Introduction reviews on the technique and the science that can result, see 
Hobbs & Dai (2017), Tiburzi (2018), Taylor (2021), and 
Verbiest et al. (2021). 

PTA data are meticulously acquired for dozens of MSPs on 


Pulsar timing arrays (PTAs) use high-precision timing of 
pulsars to attempt to directly detect and study nanohertz- 
frequency gravitational waves (GWs). The idea was first 


proposed in the 1970s (Sazhin 1978; Detweiler 1979), but it APOE days to months using the world’s largest radio 
was not until the discovery of millisecond pulsars (MSPs; telescopes. Sophisticated digital hardware synchronously 


ў m averages the faint pulsar signals at multiple radio frequencies 
Alpar et al. 1982; Backer et al. 1982) and the possibility of foma bow. bundied жерленді several ео. Айю 


calculation of times of arrival (TOAs), PTAs regularly analyze 
these data and release them to the public. 
The North American Nanohertz Observatory for Gravita- 


long-term microsecond-level timing with many MSPs that the 
idea became practical (e.g., Foster & Backer 1990). For recent 


a NASA Hubble Fellowship: Einstein Postdoctoral Fellow. tional Waves (NANOGrav; Ransom et al. 2019) was formed 
2. Physics Frontiers Center Postdoctoral Fellow. in 2007 and has previously published four data releases (see 
Deceased. 


below). Other PTA collaborations include the European PTA 
(EPTA; Chen et al. 2021), the Parkes PTA (PPTA; Reardon 
2. | et al. 2021), and the Indian РТА (InPTA; Tarafdar et al. 
и CS шде {пе terms 2022). All four collaborations compose the International РТА 


of the Creative Commons Attribution 4.0 licence. Any further 2 . 
distribution of this work must maintain attribution to the author(s) and the title (IPTA), and earlier versions of NANOGrav, EPTA, and 


of the work, journal citation and РОГ. PPTA data sets have been combined into IPTA data releases 


? NSF Astronomy and Astrophysics Postdoctoral Fellow. 
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known as IPTA DRI (Hobbs et al. 2010) апа IPTA DR2 
(Perera et al. 2019). The IPTA continues to grow, and soon 
the InPTA and MeerKAT PTA project will begin to 
contribute data to IPTA data sets (Tarafdar et al. 2022; Miles 
et al. 2023). 

This paper describes the latest release of NANOGrav data, 
the “15 yr data set,” which we have collected over more than 
15 уг using the Arecibo Observatory, the Green Bank 
Telescope (GBT), and the Very Large Array (VLA). The data 
and analysis procedures we describe are based on and extend 
those from our earlier releases, which are known as our 5 yr 
(Demorest et al. 2013, hereafter NGS), 9 yr (Arzoumanian et al. 

2015, hereafter NG9), 11 yr (Arzoumanian et al. 2018, 
hereafter NG11), and 12.5 yr “narrowband” (Alam et al. 202 1a, 
hereafter NG12) and 12.5 yr “wideband” data sets (Alam et al. 
2021b, hereafter NG12WB). We will similarly refer to the 
present 15 yr data set as NGI5. 

The data release reported on in this paper is an improvement 
on previous NG12 and NG12WB data sets in several ways. It 
includes 21 additional MSPs (for a total of 68 MSPs) and an 
additional ~2.9 yr of timing baseline for the 47 continuing 
MSPs. We have used a new Python-based timing pipeline, 
utilizing Jupyter notebooks, that relies on РТМТ” as the 
primary timing software (Luo et al. 2021). We are concurrently 
releasing both the traditional “narrowband” TOAs derived from 
many subbands of our radio observing bands and “wideband” 
TOAs derived using the methodology of Pennucci et al. (2014) 
and Pennucci (2019). We have added VLA data on six pulsars, 
one of which is a new addition to the PTA, and have included 
all of the Arecibo data available up until we were no longer 
able to use the telescope in 2020 August, 4 months before its 
tragic collapse in 2020 December. The data included in this 
release will also be included in the third data release (DR3) of 
the IPTA. 

This work is presented alongside a GW analysis searching 
for a GW background (GWB) in the 15 yr data set (Agazie 
et al. 2023c, hereafter NGl15gwb), detailed noise analysis of 
our PTA detector (Agazie et al. 2023a, hereafter NG15detchar), 
and astrophysical interpretation of the GWB results (Agazie 
et al. 2023b; Afzal et al. 2023). NANOGrav's most recent GW 
results for the NG12 data set can be found in Arzoumanian 
et al. (2020) for the stochastic GWB, Arzoumanian et al. 
(2021a) for non-Einsteinian polarization modes in the GWB, 
Arzoumanian et al. (2021b) for constraints on cosmological 
phase transitions, and Arzoumanian et al. (2023) for GWs from 
continuous wave sources. 

The plan for this paper is as follows. In Sections 2 and 3, 
we describe the observations and data reduction, respec- 
tively. In Section 4, we describe timing models fit to the 
TOAs for each pulsar, including both deterministic astro- 
physical phenomena and stochastic noise terms. In Section 5, 
we compare timing models from this data set with those from 
NGI12 and report on any newly significant astrometric, 
binary, and noise parameters. In Section 6, we acknowledge 
pulsar surveys that have discovered new MSPs included in 
this data set and present flux density measurements for these 
new pulsars at two or more radio frequencies. In Section 7, 
we summarize the work. 

The NANOGrav 15 yr data set files include narrowband 
and wideband TOAs developed in the present paper, 


73 https: / [ github.com/nanograv /PINT 
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parameterized timing models for all pulsars for each of the 
ТОА sets, configuration files used to run our timing analysis 
notebooks, and support files such as telescope clock offset 
measurements. The data set presented here is preserved on 
Zenodo at DOI:10.5281/zenodo.7967585."* Raw telescope 
data products are also available from the same website, as is 
code used to do all of the analysis described here. A living 
repository of PINT-based timing analysis software born out of 
this work can be found at https://github.com/nanograv/ 
pint pal. 


2. Observations 


The NANOGrav 15 yr data set contains data from the 100 m 
GBT, the 305 m telescope at the Arecibo Observatory 
(Arecibo) prior to its cable failure and eventual collapse, and 
the 27 25 m antennae of the VLA. The procedures we used are 
nearly identical to those in NG9, NGI1, and NGI2. As in 
NGI2, we have generated and analyzed both narrowband and 
wideband TOAs; we present both narrowband and wideband 
versions of the data set in this work. While the fundamental 
data reduction and timing procedures have not changed 
significantly from our previous releases, for this data set 
substantial effort was put into software development, with a 
major goal of improving versioning and reproducibility of 
results. Our new timing pipeline streamlines the iterative nature 
of our analyses and allows documentation of changes made to 
the software version, TOAs, timing solution, or noise model. 
The pipeline will be discussed further in Section 4. Here we 
describe the observations and data reduction used to generate 
the TOAs in this data set. 


2.1. Data Collection 


The present data set contains timing data and solutions for 68 
MSPs, 21 more than were included in NG12. A summary of the 
pulsars and their timing baselines for this data set is given in 
Figure 1. This large increase in the number of MSPs in this data 
set results from two factors: first, NANOGrav tested the 
suitability of many newly discovered MSPs and included high- 
impact PTA sources; second, NANOGrav revisited testing 
previously unsuitable (known) MSPs with new instrumentation 
for inclusion in the PTA. The MSPs with suitable timing 
precision (<1 us scatter in the daily averaged residuals over 
two to three test observations) were incorporated into the 
regular NANOGrav observing schedule prior to 2018 August, 
such that they have >2 уг timing baselines in this data set. 
Pulsars with shorter data spans are not included here, as it is 
difficult to reliably model red noise (RN) over a shorter timing 
baseline, and they do not contribute significantly to the goal of 
low-frequency GW detection. 

Data were collected using the GBT, Arecibo, and the VLA. 
A summary of the observing systems can be found in Table 1. 
The maximum timing baseline for an individual pulsar in this 
data release is 16 yr. GBT and VLA data through 2020 April 
4 are included in this data set, corresponding to the last date 
when the GUPPI back end on the GBT was used. Arecibo data 
are included through 2020 August 10, the date of the first cable 
breaking at Arecibo and the end of regular observations using 


74 All of NANOGrav's data sets are available at http: / /data.nanograv.org, 
including the 15 yr data set presented here. Raw telescope data products are 
also available from the same website. The 15 yr data set has been preserved on 
Zenodo at DOI:10.5281/zenodo.7967585. 
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Figure 1. Epochs of all TOAs in the data set. The observatory and observing frequency are indicated by color: Arecibo observations are red (327 MHZ), orange 
(430 MHz), light blue (1.4 GHz), and purple (2.1 GHz); GBT observations are green (800 MHz) and dark blue (1.4 GHz); and VLA observations are brown (1.4 GHz) 
and pink (3 GHz). The Arecibo and GBT data acquisition systems are indicated by symbols: open circles are ASP or GASP, and closed circles are PUPPI or GUPPI. 
Only a single back end (YUPPI) was used at the VLA. 


the telescope. Of the 68 MSPs in NGI5, 30 MSPs with accessible to Arecibo) with the GBT. Two pulsars, PSR J1713 
declinations 0? < 6 « 4-39? were solely observed with Arecibo, +0747 and PSR B1937+21, were observed with the GBT, 
and 31 MSPs with 6> — 45? (and outside the decl. range Arecibo, and the VLA; PSR J1600—3053, PSR 1643-1224, 
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Table 1 
Observing Frequencies and Bandwidths* 
Back Ends 
ASP/GASP PUPPI/GUPPI/YUPPI 

Telescope Frequency Usable ADM Frequency Usable ADM 
Receiver Data Span? Range? Bandwidth? Delay Data Span? Range? Bandwidth? Delay? 

(MHz) (MHz) (us) (MHz) (MHz) (us) 
Arecibo 
327 2005.0-2012.0 315-339 34 2.86 2012.2-2020.6 302—352 50 6.00 
430 2005.0–2012.3 422-442 20 1.03 2012.2–2020.6 421–445 24 1.23 
L-wide 2004.9–2012.3 1380–1444 64 0.09 2012.2–2020.6 1147-1765 600 0.91 
S-wide 2004.9–2012.6 2316–2380 64 0.02 2012.2–2020.6 1700-2404! 460 0.36 
GBT 
Revr_800 2004.6–2011.0 822—866 64 0.30 2010.2–2020.3 722—919 180 1.52 
Revrl 2 2004.6–2010.8 1386–1434 64 0.07 2010.2–2020.3 1151–1885 640 0.98 
УГА 
1.5 GHz 2015.2–2020.3 1026—2014 800 1.46 
3 GHz 2015.3–2020.0 2012—3988 1700 0.38 
Notes. 


а Table reproduced and modified from NG12WB, including adding VLA/YUPPI information. 
Dates of instrument use. Observation dates of individual pulsars vary; see Figure 1. 
* Typical values; some observations differed. Some frequencies were unusable owing to radio frequency interference. 


Approximate and representative values after excluding narrow subbands with radio frequency interference. 
© Representative dispersive delay between profiles at the extrema frequencies listed in the Frequency Range column induced by a ADM = 5 х 10 “ст ? pc, which 
is approximately the median uncertainty across all wideband DM measurements in the data set; for scale, 1 ps ~ 1 phase bin for а 2 ms pulsar with our configuration 


of пра = 2048. 
Noncontiguous usable bands at 1700-1880 MHz and 2050-2404 MHz. 


and PSR J1909—3744 were observed at both the GBT and the 
VLA; PSR J19034-0327 was observed with Arecibo and the 
VLA; and PSR J0437—4715 was observed only at the VLA. 

We adhered to a roughly monthly cadence in our observations: 
every ~3 weeks at Arecibo and every 4 weeks at the GBT and 
VLA. As described іп NG12, with the GBT and Arecibo we 
additionally performed high-cadence observations (every ~5 
days) of six of the highest timing precision pulsars for increased 
sensitivity to single sources emitting continuous GWs. There were 
some interruptions in data taking, in particular: telescope painting 
at Arecibo and azimuth track refurbishment at the GBT in 2007, 
earthquake damage at Arecibo in 2014, and a brief pause in 
observing at Arecibo following Hurricane Maria in 2017. The 
transition from test observations to incorporating new MSPs into 
our regular observing campaign at the GBT looks like a gap in 
data collection for those pulsars in early 2017. The lack of 
Arecibo data in 2018 is due to an instrumental issue described 
further in Appendix A.2. 

Each MSP was observed with at least two widely separated 
frequencies at each epoch (with a few exceptions), with the 
observing frequencies for a given pulsar chosen to maximize its 
timing precision. At Arecibo, MSPs were observed at 430 MHz 
and 1.4 GHz or at 1.4 and 2.1 GHz. Several Arecibo pulsars 
were observed at all three frequencies (Appendix B), and 
PSR J2317--1439 was originally observed at 327 and 430 MHz 
but was migrated to 430 MHz and 1.4 GHz at roughly the 
midpoint in its timing baseline. At the GBT, all MSPs were 
observed at 820 MHz and 1.4 GHz. The VLA is expected to 
provide improved sensitivity and frequency coverage in the 
2-4 GHz range, so it was used for 3 GHz observations of five 
MSPs predicted to give optimal timing precision at these 


frequencies (Lam et al. 2018b): PSR J17134-0747, PSR J1600 
—3053, PSR J1643+1224, PSR J1903+0327, and PSR J1909 
—3744. The VLA was additionally used for dual-frequency 
(1.4 and 3 GHz) observations of PSR J0437—4715, a very 
bright MSP that lies just below the decl. range accessible to the 
GBT, and for PSRJ1713+0747. The latter is observed at 
1.4 GHz with all three telescopes in order to provide a cross- 
check for possible instrumental systematics. The exceptions to 
dual-frequency observations were the high-cadence observa- 
tions at the GBT, which were carried out at only 1.4 GHz but 
with a large bandwidth. 

At Arecibo and the GBT, over the course of our experiment 
we have used two generations of data acquisition instrumenta- 
tion (i.e., pulsar “back ends"). For the first ~6 yr, the ASP and 
GASP instruments with 64 MHz bandwidth (Demorest 2007) 
were used to acquire data from Arecibo and the GBT, 
respectively, after which we transitioned to the Puerto Rican 
Ultimate Pulsar Processing Instrument (PUPPI; at Arecibo, in 
2012) and Green Bank Ultimate Pulsar Processing Instrument 
(GUPPI; at the GBT, in 2010). The PUPPI and GUPPI back 
ends processed bandwidths of up to 800 MHz (DuPlain et al. 
2008) and provided significantly improved timing compared to 
ASP and GASP. Simultaneous observations using both 
generations of instruments were used to measure offsets 
between the old and new back ends (Appendix A of NG9), 
and we continue to include these offsets in our timing models? 
as we have in previous data releases. The YUPPI back епа” 


75 These offsets are included in the ASP/GASP lines of the TOA files, after 


the “-to” flag, for pulsars with ASP or GASP data. 
76 https: / [science.nrao.edu/facilities /vla/docs /manuals /oss/performance / 
pulsar 
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was used to acquire all VLA data in this data release. The 
resulting raw data consist of folded, full-Stokes pulse profiles, 
with 2048 pulse phase bins; radio frequency resolution of 
4 MHz (ASP/GASP), ~1.5 MHz (GUPPI/PUPPI), or 1 MHz 
(ҮСРРІ); and subintegrations of 1s (with PUPPI at 1.4 and 
2.1 GHz) or 10s (with all other receiver/back-end pairings). 
Hereafter, if a description applies to all three of the GUPPI, 
PUPPI and YUPPI back ends, we will refer to them in 
shorthand as the “UPPI back ends.” 

VLA observations employed the “phased array” mode of 
telescope operation. In this mode, the voltage data streams of 
the individual antennae are coherently added, resulting in a 
single data stream with properties roughly equivalent to those 
of a single dish of equal total collecting area. This summed 
signal is then subdivided into frequency channels, coherently 
dedispersed, and folded in the same manner as typical single- 
dish pulsar data. The final results are data products (folded 
profiles) that in subsequent processing can be treated 
identically to single-dish data. In order to maximize sensitivity 
of the array sum, radio-wave phase fluctuations due to 
atmospheric variations across the array must be corrected in 
real time. During observations, phase corrections were 
measured using a bright continuum calibrator, typically within 
7-5? of the target pulsar, and applied to the data. This rephasing 
process was repeated every 10 minutes. 

For all telescopes/back ends, at the start of each pulsar 
observation, a pulsed noise diode signal was injected in order to 
later calibrate the pulsar signal amplitude during the data 
reduction steps. For GUPPI and PUPPI data, the noise signal 
amplitude was calibrated into equivalent flux density units (і.е., 
Janskys) approximately monthly using on/off observations of 
an unpolarized continuum radio source of known flux density. 
Details about the continuum sources can be found in Section 8 
of NG12. For the VLA with YUPPI, this second flux 
calibration step was not done, and the resulting data are scaled 
in units relative to the noise diode power. 


3. Data Reduction and Time-of-arrival Generation 
3.1. Data Reduction and Procedure 


We followed a standardized process to use the raw data files 
to calculate calibrated profiles and TOAs. Some steps can be 
iterative in nature; for example, initial timing must be done 
before an analysis identifying outlying TOAs. With this caveat 
in mind, we list the steps that we took to produce science-ready 
pulse profiles and to generate TOAs from those profiles: 


1. Exclusion of known bad/corrupted files or time ranges. 
From previous data sets, we were already aware of a 
number of data files that were not suitable for inclusion in 
the data set owing to various instrumental issues. In 
addition, in newer data contributing only to МС15, there 
is a ~10-month time range with unusable Arecibo data 
(see Appendix A.2). We excluded all of these files 
from NGI5. 

2. Artifact removal. This step removes artifacts from the 
interleaved analog-to-digital converter (ADC) scheme 
used by the UPPI receivers to achieve their wide 
bandwidths. The artifacts appear as negatively dispersed 
pulses and impact TOA generation if not removed. 
Details about these artifacts and their removal are given 
in Section 2.3.1 and Figure 2 of NG12. 
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3. Radio frequency interference excision and calibration. 
We performed the standard radio frequency interference 
(RFI) excision and calibration procedures detailed in 
МСОП, with the additional step of excising RFI from the 
calibration files along with the data files as described in 
Section 2.3.2 of NGI2. 

4. Time and frequency averaging. Details for Arecibo and 
GBT data are given in Section 3.1 of NG9. The UPPI data 
were  frequency-averaged  ("scrunched") into 1.6 
—32 MHz bandwidth channels, depending on the recei- 
ver. ASP/GASP data were not averaged and remained at 
the native 4 MHz resolution. The cleaned, calibrated 
profiles were time-averaged into subintegrations of up to 
30 minutes. For the few pulsars with very short binary 
periods, we time-averaged into subintegrations no longer 
than 2.596 of the orbital period. 

5. ТОА generation. This process is described in Section 3.2, 
with further details in NG9 and МО11; details for 
wideband-specific TOA generation can be found in 
NGI2WB and references therein. 

6. Outlier analysis of TOA residuals. The automated process 
by which we identified and removed outlier TOA 
residuals is described in Appendix A.4. 

7. Further, often iterative, data cleaning. Further data 
quality checks апа cleaning аге described іп 
Section 3.3. 


Steps 2-5 were performed using ће папоріре” data 
processing pipeline (Demorest 2018). This in turn uses the 
PSRCHIVE'/^ pulsar data analysis software package (Hotan 
et al. 2004) to perform most of the processing steps outlined 
here, plus PulsePortraiture for wideband template 
generation and TOA measurement (Pennucci et al. 2016). АП 
software package versions and specific data processing settings 
used in this analysis, along with the resulting TOA data (step 
5), were tracked in a git-based version control repository 
(“соадеп”). A version tag and corresponding git hash linking 
the TOA to a specific git commit are included in the final TOA 
lines, allowing a complete description of the TOA provenance 
to be reconstructed as needed. 


3.2. Time-of-arrival Generation 


TOAs are derived from folded pulse profile data by 
measuring the time shift of the observed profile relative to a 
model pulse shape for each pulsar, known as a template; this 
basic approach has been standard in pulsar timing analyses for 
decades (e.g., Taylor 1992). This data release includes both 
narrowband and wideband TOAs, as did NGI2. In the 
narrowband approach, a separate TOA is measured for each 
frequency channel, at the final time and frequency resolution 
described in Section 3.1, step 4; dispersion measure (DM) is fit 
to these TOAs along with other timing parameters. In the 
wideband approach, a single TOA and DM value are measured 
for the full receiver band, using a frequency-dependent 
template (also called a portrait) that accounts for pulse shape 
evolution (for details, see Pennucci et al. 2014). 

АП narrowband template profiles used in this data set were 
regenerated from the UPPI data, using the procedures described in 
NGS. АП profiles for a given pulsar and receiver band are aligned, 


77 https: / /github.com/demorest/nanopipe, v1.2.4. 
78 http: / [ psrchive.sf.net, v2021-06-03. 
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Table 2 
TOA Removal Flags 
Flag Reason for TOA Removal Ма М 
Initial Cuts (Loading TOAs) 
-cut simul ASP/GASP TOA taken at the same time as a PUPPI/GUPPI TOA 6883 648 
-cut orphaned Insignificant data volume (<3 epochs for a given receiver band; usually test observations) 2837 13 
-cut badrange Arecibo data affected by malfunctioning local oscillator (MJDs 57984-58447) 56,658 2036 
-cut snr Profile data used to generate TOA does not meet S/N threshold (S/N, < 8 and S/N, < 25) 255,118 1324 
-cut poorfebe Poor quality data from given front-end/back-end combination 512 5 
-cut eclipsing For pulsars showing signs of eclipses, TOAs near superior conjunction (within 10%—15% of an orbital phase) were 4551 109 
automatically excised 
-cut dmx* Ratio of maximum to minimum frequency in an observing epoch (in a single DMX bin) fj, /fmin «1.1 13,006 1086 


Initial totals 


339,565 5221 


Automated Outlier Analysis 


-cut outlier10 TOA has outlier probability р; ош > 0.1 5374 

-cut maxout Entire file removed if a significant percentage (28%) of TOAs flagged as outliers 5072 ED 
-cut epochdrop Entire epoch removed with epochalyptica (), based on an F-test p-value <10~° 6027 188 
Auto totals 16,473 188 

Manual Cuts (Human Inspection) 

-cut badtoa Remaining TOAs identified by human inspection and removed 249 1 

-cut badfile Remaining files identified by human inspection and removed 1163 170 
Manual totals 1412 171 


Removed TOAs 
Remaining TOAs 


357,450 5580 
676,465 20,290 


Notes. The flags are listed here in the order in which they were applied in our final pipeline run. A single cut flag is assigned to each TOA that is removed. 


Narrowband and wideband data sets were analyzed independently. 


а ала Г are TOA reference frequencies in the narrowband data set and are separately calculated for each TOA in the wideband data set. 


weighted by signal-to-noise ratio, and summed to create a final 
average. This average profile is then "denoised" using wavelet 
decomposition and thresholding. The process is iterated several 
times to converge on the final template. АП narrowband TOAs, 
including those from ASP and GASP data, were then generated 
using this set of templates using procedures described in more 
detail in NG9 and NG11. As in NGI2, we used an improved 
algorithm to calculate TOA uncertainties, in which we numeri- 
cally integrate the TOA probability distribution (Equation (12) in 
Appendix B of МО9) to mitigate underestimation of uncertainties 
for low-S/N TOAs. 

Wideband template portraits were created following 
NGI2WB. Again, all data for a given pulsar and receiver are 
aligned and summed, while in this case preserving (rather than 
averaging over) frequency channels. The final average portrait 
is then modeled and denoised as described in Pennucci (2019) 
to create a template that preserves profile evolution with 
frequency. For the 47 MSPs in NGI2, we reused the previous 
GUPPI and PUPPI wideband templates from that analysis. 
New templates were generated for the YUPPI data and the 21 
pulsars new to the present data set. New wideband TOAs for all 
data were then generated as described in NGI2WB. 


3.3. Cleaning the Data Set for Improved Data Quality 


In Section 2.5 of NGI2, we detailed a number of steps taken 
to improve the quality of our data set. Many of the same steps 
were taken in this data set as well, in particular the following 
data cuts (and corresponding NG12 sections): S/N cut (2.5.1), 
bad DM range cut (2.5.2), outlier residual TOA cut (2.5.3, 
except a different method was used in the present work, as 
described in Appendix A.4 below), epoch F-test cut (2.5.7), 


and orphan data cut (2.5.9). The corrupt calibration cut (NG12, 
Section 2.5.5) was not explicitly repeated in this data set, but 
the files that were cut Кот NG12 were also cut from МС15. 
Manual removal of individual TOAs or data files was done 
similarly as in NGI2, Sections 2.5.4 and 2.5.8, with some 
differences, so these cuts are described below in Appendix A.5. 

Table 2 summarizes each of the cuts used in МС15, 
including the order in which they were applied and the number 
of TOAs affected by each cut in the narrowband and wideband 
data sets. After all TOA cutting was complete, the final 
narrowband data set had 676,465 TOAs (34.646 cut) and the 
wideband data set had 20,290 TOAs (21.646 cut). 


4. Timing Analysis 


From the outset, our goal was to produce a self-consistent 
data set consisting of human-readable configuration (уат1) 
files, which produced standard timing model parameter (. par) 
files for each pulsar. Aside from the parameter values 
themselves, configuration files contain all necessary informa- 
tion to reproduce timing results, and the files are organized to 
facilitate using version control to track any changes made along 
the way. All timing fits used the JPL DE440 solar system 
ephemeris (Park et al. 2021) and the TT(BIPM2019) timescale. 

As in NG12, we used standardized Jupyter notebooks 
(Kluyver et al. 2016) to automate our timing procedures and 
built new software tools to improve transparency, code 
readability, and reproducibility of our results. Our analysis 
pipeline is now entirely PINT based (Luo et al. 2021)— 
previous releases used TEMPO”? (Nice et al. 2015)—and 


7? https: / / github.com/nanograv Летро 
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flexible enough to accommodate working with both narrow- 
band and wideband TOAs. A frozen version of code used to 
produce Ма15 results, as well as intermediate data products, 
can be found at data.nanograv.org (see Section 4.5 for more 
details). An open-source version of our timing analysis 
package, PINT Pal, is also available.?? 


4.1. Timing Models and Parameters 


Pulsar timing requires determination of the TOAs and 
comparison to a model composed of many different physical 
effects. Our timing parameter sets fall into six categories, with 
the numbers of parameters in each category fit per pulsar shown 
in Appendix B. 


1. Spin: We fit three parameters describing the rotational 
phase, frequency, and frequency derivative of each 
pulsar. In one case (long-period binary PSR J1024 
—0719), we also fit the second frequency derivative. 

2. Astrometry: For each pulsar we fit five parameters 
describing the two-dimensional position on the sky in 
ecliptic coordinates, the two-dimensional proper motion, 
and parallax. As in our previous works, we fit all five of 
these regardless of the measurement significance. 

3. Binary: For 50 of our observed systems we fit binary 
parameters describing the orbit with a companion star. 
With the exception of the long-period binary PSR J1024 
—0719, we fit, at a minimum, the five Keplerian 
parameters fully characterizing the orbit. Several different 
models were used depending on the orbital character- 
istics; see Section 4.1.1. 

4. Dispersion measure: 'Time variations in the DM require 
that we fit the TOAs with a dispersive delay proportional 
to v ? and the unknown DM at that epoch. This is 
discussed in greater depth in Section 4.2. 

5. Frequency dependence: Additional time-independent but 
frequency-dependent delays are fit for on a per-pulsar 
basis (“FD” parameters; see NG9), with the number 
included determined via our F-test procedure discussed 
below in Section 4.1.2. In the narrowband data, we take 
these to describe time offsets due to differences in the 
observed pulse shape at a given frequency compared to 
the template shape used in timing. In the wideband data, 
we expect the pulse portrait to encapsulate the changing 
frequency dependence of the profiles. However, we still 
find several significant FD parameters. We fit for these 
values in the same manner regardless, and the signifi- 
cance of these will be explored more in future work. 
Though the physical mechanism has not been definitively 
determined, other time-varying chromatic processes 
could be picked up as time-independent FD parameters. 

6. Jump: We fit for "jumps" to account for unknown phase 
offsets between data observed by different receivers and / 
or telescopes. 


A complete set of timing residuals and DM variations for 
each pulsar in our data set can be found in Appendix C. 


4.1.1. Binary Models 


We used one of five binary models, depending on the orbital 
characteristics of the pulsar in question. Low-eccentricity orbits 


20 https: //github.com/nanograv /pint_pal 
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(see Section 5.1.3) were fit using the ELL1 (Lange et al. 2001) 
or ЕЛІН (Freire & Мех 2010) model, depending on whether 
or not Shapiro delay was marginally detected; orbits with 
higher eccentricity were modeled with the BT (Blandford & 
Teukolsky 1976) or DD (Damour & Deruelle 1985, 1986; 
Damour & Taylor 1992) model; and for PSR J17134-0747, in 
which we measure the physical orientation of the orbit, we used 
the DDK model (Kopeikin 1996). Pedagogical descriptions of 
these binary models can be found in Lorimer & Kramer (2012). 
In all cases, at a minimum we included and fit five Keplerian 
parameters in the binary model: the orbital period P; or orbital 
frequency fp, the projected semimajor axis x; and either the 
eccentricity e, longitude of periastron w, and epoch of 
periastron passage То in the DD or DDK models, or two 
Laplace-Lagrange parameters (еі, 62) and the epoch of the 
ascending node Т in the ЕПЛ ог ELLIH models. 

The ELL1 model describes the orbit with the Гарјасе– 
Lagrange eccentricity parameterization, with а = esinw and 
еҙ = ecosw (for a more complete description of this parameter- 
ization, see Lange et al. 2001). This parameterization is needed 
because for a nearly circular orbit one cannot reliably define the 
time and location of the periastron. If the Shapiro delay was 
marginally detected via the F-statistic as described below, we 
instead modeled low-eccentricity orbits with the ELLIH model, 
which is simply the ELL1 model with the addition of the h3 and 
h4 parameters from the orthometric Shapiro delay parameteriza- 
tion of Freire & Wex (2010). For a significant detection of 
Shapiro delay, the companion mass rn, and orbital inclination sin i 
parameters were instead included in the ELL1 model. 

The BT and DD models directly measure e, x, and То and аге 
more accurate than ће ELL1 model for orbits with a large enough 
eccentricity to measure the longitude and epoch of periastron. The 
BT model is Newtonian, while the DD model is a theory- 
independent relativistic model that allows for the inclusion of 
Shapiro delay parameters, companion mass m, and inclination 
angle via sin i. In our data set, the use of ELL1 and DD is split 
nearly evenly between pulsars. For one of our most precisely 
timed pulsars, PSR J17134-0747, we found in NG12 that we were 
sensitive to the annual-orbital parallax (AOP; Kopeikin 1995), 
which allows the true orientation of the orbit to be measured. 
Thus, we transitioned from the DD to DDK model for this pulsar. 
The DDK model incorporates the orbital inclination and the 
longitude of the ascending node of the orbit, Ол (Kopeikin 1996). 
Appendix D describes our use of the DDK model in detail. 

For the majority of orbital models, we used the orbital period 
Pp, but for four short orbital period (P, < 0.5 days) pulsars we 
instead included the orbital frequency f; in the model. Using the 
orbital frequency instead of period allows us to include multiple 
orbital frequency derivatives to better describe the orbit, rather 
than being restricted to only the first period derivative if using РЬ; 
it also allows us to easily test for additional derivatives using the 
F-statistic. All four of these pulsars are in “black widow" systems 
(Swihart et al. 2022). Three of them are noneclipsing (PSR J0023 
+0923, PSR J0636+5128, and PSR J22144-3000), while 
PSR J1705— 1903 eclipses for ~30% of each 4.4hr orbit. All 
four of these pulsars have low-eccentricity orbits and were 
modeled with the ЕПЛ model. For PSR J2214+3000, only }, 
(“ЕВ0” in the timing model) was needed, with no derivatives; 
J06364-5128 required three derivatives (FBO through FB3); and 
10023--0923 and J1705—1903 required five derivatives (ЕВО 
through FB5). 
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In all cases, we tested for the presence of the Shapiro delay 
using the F-statistic, as described in Section 4.1.2 below. If m, 
and sin i were significantly detected, we added these parameters 
to the timing model and refit. 


4.1.2. Parameter Inclusion and Significance Testing Criteria 


As in NG12 and our prior data releases, an F-statistic was 
used to determine what additional parameters should be 
included in the timing model beyond those describing pulsar 
spin, astrometry, RN and white noise (WN), and Keplerian 
binary motion. This test is valid in the case that one model is a 
subset of another; therefore, these tests were performed by 
fitting a nominal model, recording the resulting x^, then adding 
or subtracting the parameter of interest, and refitting. The F- 
statistic is 


(о x2)/@o = п) 

х2/п 
where х? and по are the best-fit chi-squared value and number of 
degrees of freedom, respectively, for the nominal (superset) 
model. Parameters were included if they induced a р < 0.0027 
(~За) change in the fit. As described in NG12, this test was 
applied to secularly evolving binary parameters (В), c», and x), 
Shapiro delay parameters, frequency-dependent (FD; see NG9) 
parameters, and higher-order orbital frequency derivatives (ЕВл, 
where n=O corresponds to the orbital frequency and n = (1, 
2,...} correspond to frequency derivatives) where applicable. For 
parameters with components of increasing complexity (like 
FD1-5, where a model cannot include noncontiguous combina- 
tions of parameters), the test was applied iteratively. All FB and 
FD terms below the highest-order term that induced a significant 
change in the fit were included, even if a subset of the lower-order 
parameters did not produce a significant change; for example, if 
FD3 is significant but FD2 is not, then FD1, FD2, and FD3 are 
included in the model. The summary plots generated with 
PINT Pal provide the user with suggestions for parameters to 
add or exclude in the model. 


F 


(1) 


4.2. Dispersion Measure Variations 


Changes in the Earth-pulsar line of sight require that we 
account for variations in the DM on short (i.e., per-epoch) 
timescales (Jones et al. 2017). We measure the pulsar signals 
over a wide range of radio frequencies to estimate the time- 
varying DM for both narrowband and wideband data sets (see 
Figure 1). We used the DMX model, a piecewise constant 
function, to describe these variations, with each DMX 
parameter describing the offset from a nominal fixed value. 
Modeling of these variations requires up to hundreds of 
additional parameters (see Appendix B), and for most of our 
pulsars the expected variations from interstellar turbulence 
alone are significant enough between epochs that we cannot 
reduce the number of parameters (Jones et al. 2017). 

At Arecibo and the VLA, observations with different 
receivers were consecutive, whereas with the GBT one or 
more days passed between observations for scheduling 
efficiencies.*! We assumed a constant value for the DM over 
a window of time 0.5 days long for any pulsars with Arecibo 


B Switches between prime-focus receivers (i.e., “Ксуг 800”) and Gregorian- 
focus receivers (i.e., *Rcvrl 2"), or vice versa, at the GBT take 10—15 minutes. 
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observations and longer windows for pulsars solely observed 
with the GBT, with a 15-day range for GASP data and 6.5-day 
range for GUPPI/YUPPI data. If the maximum-to-minimum 
frequency ratio of the TOAs did not satisfy va, / Ини > 1.1, we 
excised the data as described above and did not fit а DMX 
parameter. 

Besides DM variations due to the turbulent interstellar 
medium, another visually apparent contribution is from the 
Earth- pulsar line of sight intersecting different sections of the 
solar wind. For a few specific pulsars close to the ecliptic plane, 
the change in DM can be significant even within the time range 
of a single window. We use the following criterion to find such 
windows. As in NG11 and NGI2, for this purpose we used а 
toy model of a spherically symmetric and static solar wind 
electron density given by n,(r) = no(r/ro) 2, where г is the 
distance from the Sun and ny = 5 ст > is the fiducial electron 
density at a distance го = 1 au (e.g., Splaver et al. 2005). If the 
projected DM variation due to the change in the line of sight 
through the solar wind would have induced a timing variation 
of greater than 100 ns, the DMX time ranges were divided into 
0.5-day windows, with data still subject to the frequency cut 
described above. Note that this toy solar wind model is used 
only for this test while setting DMX window sizes; it is not 
included in the pulsar timing models themselves. Thus, 
reported DMX values in the timing models incorporate 
dispersion by both the solar wind and the interstellar medium. 

In addition to being our most-observed pulsar, PSR J1713 
+0747 experienced rapid chromatic variations over a period of 
several months in 2016 (e.g., Lam et al. 2018a; Chen et al. 
2021; Goncharov et al. 2021) such that it was necessary to split 
the DMX windows to avoid additional excess timing noise 
from the event as in №012. Following the conclusion of the 
current data set, PSR J1713+0747 experienced a profile 
change, first reported in Xu et al. (2021), confirmed by Meyers 
& Chime/Pulsar Collaboration (2021), and further discussed 
by Singha et al. (2021) and Jennings et al. (2022). As this event 
occurred outside of the range of data reported here, we do not 
need to correct for it in this release. 


4.3. Noise Modeling 


Pulsar timing uses a phenomenological noise model for the 
data that is separated into two broad components distinguished 
by the timescale of the correlations between the TOAs. The 
WN model, i.e., noise that is uncorrelated in time and hence 
spectrally flat, inflates the values of the TOA uncertainties, 
Озум, derived from the pulse template matching process 
(Taylor 1992). The template matching process assumes that 
the data are a scaled and shifted version of the pulse profile; 
however, other sources of noise can be best modeled as WN 
but do not come from the template matching process (Lommen 
& Demorest 2013). 

Our noise modeling paradigm is discussed in detail in 
NGl5detchar and briefly summarized here. Three WN 
parameters—EFAC (F), EQUAD (О), and ECORR (7)— 
are used to adjust TOA uncertainties in order to accurately 
represent the WN present in the data. EFAC scales TOAs 
linearly, accounting for uncertainty induced by template 
matching errors or template mismatches. For well-characterized 
systems, EFAC tends to 1.0. EQUAD adds an additional WN 
in quadrature, ensuring a minimum error size. ECORR also 
adds in quadrature and accounts for uncertainty that is 
correlated among frequencies within an observation, most 
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importantly pulse jitter, though other mechanisms can add as 
well (Lam et al. 2017). Various differences between pulsar 
timing back ends and radio observatory receivers (e.g., the 
frequency ranges covered) make it necessary to give different 
values of these WN parameters to each receiver/back-end 
combination. These WN terms come together with receiver/ 
back-end combination (re/be) dependence as 


Су = F(re/be) [ogni + Q'(re/be)] бу + 7 *(re/be) Му, 
(2) 


where the i, j denote TOA indices across all observing epochs, 
ду is the Kronecker delta, and we omit the dependence on 
receiver and back end, re/be, from here on for simplicity. 
While EFAC and EQUAD only add to the diagonal of C, 
where С; are the elements of the covariance matrix, the 
ECORR terms are block diagonal for single observing epochs. 
ECORR is modeled using a block diagonal matrix, М, with 
values of 1 for TOAs from the same epoch and zeros for all 
other entries. 

Pulsar timing data often show evidence of correlations across 
longer timescales than can easily be modeled with ECORR 
terms. These correlations are instead modeled as a single 
stationary Gaussian process with a power-law spectrum 


ШЕТІ 
мр = ла]. 
lyr 


where f 15 the Fourier frequency, Yea is the spectral index, 2 
and Ага is the noise amplitude at a reference frequency of 
1 yr !, Long-timescale correlated data are characterized as RN 
in which 4444 « 0. ECORR is not separately modeled in the 
wideband data set, as it is completely covariant with EQUAD. 

The timing model and noise analysis are performed 
iteratively over each pulsar data set, first fitting a preliminary 
pulsar ephemeris and then doing a Bayesian noise analysis with 
ENTERPRISE? using a linearized timing model (see 
Section 4.4 below for details about the linearized version of 
the timing model). We then refit the timing model with the 
noise parameters obtained from ENTERPRISE. This process is 
repeated until the timing and noise parameters stabilize and 
noise posteriors look nearly identical from one run to the next. 
Slight changes from NGI2 in the implementation of the 
Bayesian noise analysis include increasing the number of 
samples obtained and matching the prior for the spectral index 
to those used in the GW detection pipeline. 

The WN parameters EFAC, EQUAD, and ECORR are 
measured for each pulsar/back-end/receiver combination. A 
small fraction of these WN parameters changed by >30 
between NGI2 and NGI5 (28/159, 17/159, and 5/159 for 
EFAC, EQUAD, and ECORR, respectively). Even so, these 
changes are not unexpected as the length of the data set 
increases, especially for pulsars that were newly added 
to NG12. 

Twenty-three pulsars in NGI5 were found to have 
significant levels of RN. Ten of these measurements are newly 
significant in NG15, while 13 of the 14 sources with significant 
RN in NGI2 continued to favor its inclusion in their timing 


(3) 


82 Note that in this paper Yreq is the true spectral index of the RN spectrum and 
could be positive or negative. In other papers this equation is often specified to 
have a negative spectral index such that the value of “уга is positive. 

83 https: //github.com/nanograv /enterprise 
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models. Only PSR J23174-1439, which favored an Аа an 
order of magnitude lower than the next-lowest Аа in NG12, as 
well as a steep “га comparable to that of PSR J00304-0451, 
had RN parameters removed from its timing model in NGI5. 
No еа OF Аа values differed between the data sets by more 
than a factor of ~3а for pulsars with significant RN in both 
NG12 and NGI25. 


4.4. Linearity Testing for Pulsar Timing Models 


A pulsar timing model is a deterministic function #4 that 
takes parameter values p and predicts observed pulse phases at 
the time of each TOA 199 (in time units) such that the observed 
TOAs are the deterministic values plus noise и, i.e., 
[9s = уде (раде) ^ п. We are able to simplify our GW 
detection calculations by using a linear approximation to this 
function (Ellis et al. 2013): 


Ot (Prue) 


др + п, 
Обр Ши 


r (Poest) = p= 1*9 (pies) ~ 
ép =0 


(4) 


where Ppest is our vector of best-fit parameters from our 
generalized least-squares fit that form our residuals r, and 
бр = Prue — Poest 18 therefore small. The term with the partial 
derivative on the right-hand side is the design matrix. One 
generally expects this approximation to be good when the data 
constrain the parameters tightly, relative to the level of 
nonlinearity intrinsic to each parameter; for example, some 
are entirely linear. Since this level varies significantly between 
different parameters in the model, we developed tests to verify 
that the approximation of linearity is adequate. 

In our reported models, the parameter uncertainties and the 
parameter correlations define an ellipsoid. If the linear model is 
adequate, then the plausible pulsar parameter values are 
constrained by the data to lie within or near this ellipsoid. 
Thus, this ellipsoid is the region where we require the timing 
model to be well approximated by the linearized model. While 
we generally use analytical marginalization with flat priors to 
remove the (linearized) timing model parameters from the GW 
detection process, we expect that most samples from 
preliminary detection Markov chains will fall within or near 
this error ellipsoid (Kaiser et al., in preparation). Thus, we try 
to test the degree to which the full timing model differs from its 
linearized approximation over this error ellipsoid. 

One direct approach is simply to shift the value of each 
parameter, one at a time, by 1c from the best-fit рь. For each 
parameter, this produces a trial parameter set риа, from which 
we Шеп compute Ше difference in residuals 
T (Puia) = (риа) (in time units). This provides а measure- 
ment of the deviation from linearity for this parameter at the 
time of each TOA. In the case of highly covariant parameters 
this approach may result in p,;4 falling outside the full 
(multidimensional) error ellipse. Therefore, this test can be 
considered a conservative upper limit on the level of 
nonlinearity; the true level may be lower. For each pulsar 
and parameter we checked whether the rms amplitude of this 
deviation was less than 100ns. We found only a few 
parameters on a few pulsars that exceeded this threshold. The 
most common example was the time (То) and longitude (ш) of 
periastron. These two parameters form part of an essentially 
polar representation of the orbital shape, and when the orbit is 
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nearly circular, this coordinate system introduces nonlinea- 
rities; for the pulsars where this was a problem the nonlinearity 
was typically on the scale of a few microseconds. A second 
cause for nonlinearity is the parameter sini, which is 
constrained by the Shapiro delay. The shape of the Shapiro 
delay as a function of orbital phase is quite nonlinear in the 
parameter sini, and so for pulsars where the Shapiro delay is 
poorly constrained this nonlinearity introduces a few-micro- 
second discrepancy. 

For the pulsars and parameters where nonlinearities appear 
in the above test, we carried out a more stringent test. Since the 
detection tools analytically marginalize the timing parameters, 
any discrepancy that can be modeled as a linear combination of 
the timing model derivatives disappears. We therefore 
evaluated each of the problem cases, removing all such 
components. The discrepancy that remains is not absorbed 
during Bayesian fitting. This process resolved nearly all of the 
above nonlinearities, with the following exceptions: the 
parameter sini for the pulsars J1630+3734, J1811—2405, 
J1946--3417, J20174-0603, and J2302+4442, and the para- 
meters sin i, Рь, ш, ш, and То for the pulsar J1630--3734. For 
the parameter sini, in each case the Зо confidence interval 
extends above the physical limit of 1, so nonlinearities must be 
expected. They take the form of a sharp peak around superior 
conjunction and are mostly about 1 jus. PSR J1630+3734 is a 
pulsar with only 3 yr of data, and because of the orbital 
coordinate system, its P, is highly correlated with w, along with 
the nonlinearities we observe in other pulsars with nearly 
circular orbits. The discrepancies in this case are of the order 
10 us. 

The conclusion of this analysis is that the vast majority of 
pulsar parameters in our data set can be very well approximated 
by a linearized model; therefore, applying this approximation 
when marginalizing over the timing model as part of GW analyses 
is not a significant source of error. For the small number of 
parameters that may show significant nonlinearity, the effect 
occurs at the pulsar's orbital period. These are generally days to 
weeks and so are unlikely to interfere with the GW signals of 
interest on timescales of years. The main result of nonlinearity is 
that parameter uncertainties estimated from the linearized model 
are potentially unreliable. For these parameters additional analysis 
is required to determine posterior probability distributions; refer to 
Section 5.1.4 for discussion of Shapiro delay measurements and 
pulsar masses in our results. 


4.5. Data Products 


Our complete catalog of observations for 68 MSPs is 
available alongside this publication at https:/ /data.nanograv. 
org. We also provide metadata and other useful intermediate 
files generated from pipeline runs there, including the 
following: 


. configuration (*. yam1) files; 

. standardized pulsar parameter (* . par) files; 

. initial TOA (*. tim) files; 

. TOA files with cut flags applied ('excise.tim) 
. noise modeling chains and parameter files; 

. DMX time series ('dmxparse".out). 


холоо мыо н 


For ease of use, timing models will be made available in both 
Tempo2** (Hobbs & Edwards 2012) and PINT formats. In 


84 https: //bitbucket.org /psrsoft /tempo2 
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addition, the correlation matrices for fit parameters of all 
pulsars are provided in the following formats: 


1. human-readable list of pairwise correlations (*. txt); 

2. machine-readable NumPy compressed correlation 
matrices (*.npz); 

3. standardized HDF5 compressed correlation matrices 
(^. hd£5). 


Finally, we have containerized the production environments 
utilized in the production of NG15. The following options will 
be released alongside the NG15 data set: 


1. An interactive Docker (Merkel 2014) container, which is 
built off of the Official Jupyter images to ensure 
compatibility with the JupyterHub software stack (Gran- 
ger & Pérez 2021). 

2. A Singularity (Kurtzer et al. 2017) container, which is 
optimized for HPC workloads. 

3. An Anaconda Environment, consisting of Conda-Forge 
(conda-forge Community 2015) packages tested against 
Python 3.10. 


5. Timing Model Comparisons and Newly Measured 
Parameters 


5.1. Comparison of NG12 and NGI5 Timing Models 


The careful analysis of changes to timing model parameters 
can provide a wealth of secondary science results and provide a 
measure of confidence in addition to the data set. Here we 
present a basic summary of changes between the 12.5 and 15 yr 
data releases, focusing specifically on improvements in 
astrometric and binary parameter measurements, changes in 
WN and RN parameters, and discrepancies between timing 
solutions obtained with Tempo2 and PINT. 


5.1.1. Astrometric Parameter Comparison 


We begin by comparing the astrometric parameters, 
specifically the parallax and proper-motion values. We are 
particularly interested in NGI5 parameters that have changed 
by >30 2.5 relative to NG12, where 6125 is the 1c uncertainty 
on the NG12 parameter (i.e., a conservative estimate of a 3012.5 
discrepancy). We find three pulsars with 73015 5 changes in the 
ecliptic longitude component of proper motion и; (PSR J0636 
+5128, 3.50125; PSR J1909—3744, 3.70155; and PSR 11946 
+3417, 7.004,55) and one pulsar with a 30,55 change in 
parallax (PSR J00234-0923, 3.4125). It is known that varia- 
tions in the RN parameters, orbital parameters, and other 
parameters that require long timescales for measurement have 
covariances with astrometric parameters, such that a few 
pulsars will have 23е changes in some of these parameters. 
Such covariances can therefore explain the changes in 
astrometric parameters between МО12 and NGI5. For 
example, PSR J1946+3417 has RN in NGIS5 but not in 
NGI2, explaining the change in proper motion, and the 
measured parallax of PSR J00234-0923 changes significantly 
with the addition or removal of orbital frequency derivatives in 
its timing model. 

In all cases, the uncertainty on the 15 yr parameter is smaller 
than that of the 12.5 yr parameter, as is expected for а 
parameter measured with a longer data set. Figure 2 shows the 
12.5 and 15 yr parallax values for all pulsars that were included 
in both data sets. 
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Parallax with l-o uncertainties 


J0023--0923 
J0030--0451 
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J1024—0719 
J1125+7819 
J1453+1902 
J1455—3330 
J1600—3053 
J1614— 2230 
J1640--2224 
J1643—1224 
J1713+0747 
J1738+0333 
J17414-1351 
J1744—1134 
J1747—4036 
J1832—0836 
J1853+1303 

В1855--09 
J1903+0327 
J1909—3744 
J1910+1256 
J1911+1347 
J1918—0642 
J19234-2515 

B1937+21 
J1944--0907 
J1946+3417 

B19534-29 
J2010—1323 
J20174-0603 
J2033+1734 
J2043+1711 
J2145—0750 
J2214+3000 
J22294-2643 
J2234+0611 
J22344-0944 
J23024-4442 
J2317+1439 
J2322+2057 


Parallax (mas) 


Figure 2. Measurements and uncertainties of parallax from the 12.5 yr and the current 15 yr data sets, showing the parallax values to be consistent across data sets. 
Orange crosses denote the NG12 measurement, while blue circles denote the NG15 measurement. Horizontal lines of the corresponding color mark the extent of the 
lo measurement uncertainty. The timing algorithm allows both positive and negative parallaxes, even though only positive values are physically meaningful. As 
expected, none of the negative parallax values are significant within their uncertainties. 


5.1.2. Comparison of Narrowband and Wideband Timing Models interpretation-related results. Its primary purpose is to aid in 
As in NG12 and NGI2WB, we have generated both a refining the techniques used in wideband TOA generation and 
narrowband data wet and а wideband data set as part of the in the curation of wideband timing solutions in order to prepare 


NG15 process. However, it is critical to note that the wideband for the eventual inclusion of data from wide-bandwidth 
data set is not being used to derive any GWB or astrophysical- receivers and cope with prohibitively large narrowband data 
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volumes in future releases. In this work, the wideband data set 
is only used to derive preliminary pulsar masses for sources 
with significantly measured Shapiro delay parameters (see 
Section 5.1.3). 

A thorough comparison of all measured parameters common 
between the narrowband and wideband timing solutions for 
each source was conducted, as in NG12WB. In that work, few 
discrepancies larger than 20125 were found. However, the 
longer NGI5 wideband data set has shown a handful of 
discrepancies that are currently under investigation. In several 
cases, the length of either the wideband or the narrowband data 
set is sufficiently long (due to automatic outlier excision in the 
other data set) to significantly measure RN. In others, either the 
data set length or presence of RN caused a parameter (usually 
orbital) to be included in one of the data sets, but not both. 
Covariances between (usually orbital) parameters as discussed 
in Section 4.4 may also worsen these discrepancies. For 
example, a barely significant measurement of the relativistic 
Shapiro delay would require the addition of the two post- 
Keplerian parameters that describe the effect. However, if the 
constraint is not strong, nonlinearities in binary orbital 
parameters may exacerbate discrepancies in these parameter 
measurements. In general, these covariances might lead to 
larger possible discrepancies in narrowband and wideband 
parameter constraints. The overall impact of these discrepan- 
cies is minimal, as the only analysis referencing the wideband 
data set is the preliminary measurement of pulsar masses. 
Further refinement of this procedure will precede the publica- 
tion of the next NANOGrav data release. 


5.1.3. Binary Model Comparison 


The high cadence and long timing baseline of NANOGrav's 
data set enable the detailed study of a large number of MSP 
binary systems. NG9 and its accompanying publication 
Fonseca et al. (2016), as well as NG11 and NGI2, have 
included analyses of binary sources timed as part of the 
NANOGrav PTA. In this work, we present a basic overview of 
newly measured binary parameters, changes to timing models, 
and new constraints on pulsar masses from the measurement of 
post-Keplerian Shapiro delay orbital parameters. A manuscript 
describing our more in-depth analysis of the 15 yr binary 
sources is in preparation. 

Of the 68 MSPs included in the current data set, 50 are in 
binary systems. Because these sources are selected for their 
high timing precision and stability, the vast majority are in 
near-circular orbits with white dwarf companions. Of the 21 
new MSPs added since МО12, 18 аге binaries. As more 
observations are added to the data set, sensitivity to secular 
binary parameters (i.e., X, ш, and P,) improves. Additionally, 
parameters such as the relativistic Shapiro delay can be more 
precisely measured with better orbital phase coverage. We see a 
number of improvements and changes to binary parameter 
measurements resulting from the addition of —3 yr of 
NANOGrav observations. 

Consistency between the narrowband and wideband data sets 
is discussed in Section 5.1.1. Any discussion regarding the 
addition and removal of parameters, as well as their values, is 
based on the narrowband NG12 and 15 yr data sets. However, 
because the wideband data set consists of significantly fewer 
TOAs (the NG12WB data volume is ^33 times smaller than 
that of NG12), mass determinations based on post-Keplerian 
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Table 3 
Binary Models and Binary Model Changes for MSPs Common to NG12 

and NGI5 

Source Binary Model New Parameters Removed Parameters 

Ј0023+0923 _ ELLI ЕВ5 © 

J0613-0200 X ELLI1 P, x 

J0636+5128 | ЕПЛ FB2, FB3 

J0740--6620 ЕШЛІ P, 

J1012+5307 — ЕПЛ x 

1125-7819 ELLI x 

11455-3330 рр Ue 

Ј1600–3053 рр Ф 

11614-2230 ЕШ 1 Us 

J1640+2224 DD é 

J1643-1224 DD é 

1713-0747 DDK 

1738-0333 _ ЕПЛ Ue 

11741-351 ELLI РЬ 

1853-1303  ELLIH to DD 

В1855--09 ЕПЛ 

Ј1903+0327 рр 

11909-3744 ELLI 

J1910+1256 DD Ue 

Ј1918–0642 ЕПА to DD % 

11946+3417 DD 

В1953--29 рр 

Ј2017+0603  ELLI 

J2033+1734 DD ee 

J2043+1711 ЕМІ РЬ 

J2145-0750 ELLIH 

Ј2214+3000  ELLI 

J2229+2643 DD Е 

J2234+0611 рр РЬ 

J2234+0944 ELLI 

J2302+4442 рр 

J2317+1439  ELLIH 


parameter measurements are based on the wideband data set to 
reduce computational burden. 

Discrepant binary parameters: Several NANOGrav MSPs 
have shown significant changes (>30) in measured Keplerian 
orbital parameters since МС12, as determined by the absolute 
value of the difference in parameter value between NGI2 and 
now, divided by the NG12 uncertainty. However, these >30 
discrepancies can generally be explained by changes to the 
pulsar's orbital or RN model. 

For example, the significance of P, for PSR J1600—3053 
decreased since NG12 primarily owing to the addition of RN, a 
new measurement of w, and an added frequency-dependent 
parameter (FD3). Relative to the recent, larger uncertainty, the 
value of P, increased by 4.30. Newly measured RN in 
PSR J1614—2230 (one of the highest-mass neutron stars 
known; see Demorest et al. 2010) was accompanied by a 
3.8c increase in Рь. Such changes are plausible as the РТА 
timing baseline increases and sensitivity to long-term orbital 
variations improves. Because secular evolution of parameters 
can happen on timescales similar to those of intrinsic RN, these 
two measurements can be covariant (see NG12). 

Ten of the 50 binary sources show newly measured values 
for either P,, x, or ё (see Table 3). In the case of PSR J0613 
—0200, х was removed and Б, was added. 

Changes to binary models: Because the ELL1 binary model 
is only valid for sufficiently circular orbits, a test 
(a sin(i)e?/ C «€ Og IN / NA; see Appendix A1 of Lange 
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et al. 2001) was imposed on all ELLI binaries to ensure its 
validity. ЕШ 1 is no longer valid for PSR J1918—0642, so this 
source is now parameterized by the DD binary model. 

The sources PSRJ1853--1303 and PSRJ1918—0642 
required changes to their binary models as a result of the 
increased length of the data set; specifically, we have chosen to 
use the DD binary model instead of ELL1H and ELL1 for PSR 
J1853+1303 and PSR J1918— 0642, respectively (see Table 3). 
PSR J1853--1303 was discussed in detail in NG12 because the 
orthometric Shapiro delay parameters Ла and h4 were newly 
constrained. We now significantly measure the traditional 
Shapiro delay parameters m, and sini as determined by an F- 
test comparison; therefore, the DD parameterization is merited. 
Although these parameters are significantly measured, we do 
not yet meaningfully constrain the pulsar mass (my = Bau 
Мо; see Section 5.1.4). Continued observations will result in 
improved orbital coverage, potentially enabling a more precise 
my measurement in future data releases. 


5.1.4. Pulsar Masses from the Relativistic Shapiro Delay 


Because an extensive analysis of each of the 50 binary MSPs 
included here is beyond the scope of this data release, a more 
in-depth manuscript detailing our binary analysis similar to 
NG9 and Fonseca et al. (2016) is in preparation. Here we 
present pulsar mass (m,) measurements obtained from 
measurement of the relativistic Shapiro delay in cases where 
the companion mass (те) and orbital inclination angle (i) are 
significantly constrained. When these two parameters are 
combined with the Keplerian mass function and the extremely 
well-determined x and Р, orbital parameters, the pulsar's mass 
and companion mass can be measured independently. In the 
rare case that other post-Keplerian parameters can also be 
constrained, one's ability to precisely measure m, is improved; 
however, our basic analysis does not take this additional 
information into account. 

Posterior probability distribution functions (pdf's) for pulsar 
masses (see Figure 3) were derived from grid-based iteration 
over т„ and sin(i) using PINT. For each combination of m, and 
sini, a x^ fit was performed without holding other model 
parameters (except the measured WN and RN values) fixed. 
Equations (4) and (5) and associated text in Fonseca et al. 
(2016) explain the Bayesian translation (with uniform priors) 
from this x grid to marginalized posterior pdf's in т. and sin i 
in greater detail. While Fonseca et al. (2016, 2021, 
hereafter, FCP21) perform the x model fits using Tempo2, 
the PINT-based results presented here were cross-checked with 
Tempo2-derived results. For assumed white dwarf companions 
with poorly constrained masses, grids spanned the 0-2 Mo 
range. АП inclination angles (sampled uniformly in cos і rather 
than sini to represent a random distribution of orbital 
orientations) were searched if that parameter was not already 
well constrained. Unless otherwise noted, reported uncertain- 
ties correspond to 68.3% confidence intervals, and grids were 
200 by 200 samples in size. 

Two of the 18 binary MSPs added to the NANOGrav data 
set since NG12, PSR J1630+3734 and PSR J1811—2405, 
show (according to the F-test) significantly—if not precisely— 
measured Shapiro delay (see Table 4). We measure the mass of 
PSR J1630+3734 to be 6.7439 Мо. While this mass is 
intriguingly high, this source has only 3 yr of timing data 
and suffers from nonlinearities in multiple of its measured 
binary parameters (see Section 4.4). PSR J1811—2405 suffers 
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12.5-year vs. 15-year Measurements of mp 
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Ј1946+3: 


Ј2017 


J2302+4 


Pulsar Mass mp (Mo) 


Figure 3. Posterior probability distributions for pulsar mass m, for each NG15 
binary with well-constrained traditional Shapiro delay parameters. Solid lines 
within the pdf curves indicate the median mass value, while the two dotted 
lines represent the 68.3% confidence interval. Orange curves were derived 
using NG12 measurements; blue curves are from the updated NGI5 data set. 
Of the 14 Shapiro delay measurements common to NG12 and NGI5, seven of 
the mean masses decreased, while seven either increased or did not appreciably 
change. 


from the same limitations and nonlinear tendencies as 
PSR J1630+3734. A lengthened timing baseline will help 
resolve some of these covariances and improve our constraint 
of m,. 

А number of factors determine the precision of Shapiro delay 
measurements. Among these are the density of observational 
coverage around superior conjunction, a pulsar’s timing 
precision and noise characteristics, and its orbital geometry. 
A number of m, measurements have improved between NG12 
and the current data release, due in part to a longer timing 
baseline and improved orbital coverage (see Figure 3 and 
Table 4). The 68.3% m, confidence intervals for PSR J0740 
+6620, PSR J1600—3053, and PSR J17414-1351, improved 
by a factor of >2.5 between NG12 and the present work. 

PSR 10740--6620 is of particular interest to those wishing to 
constrain the poorly understood dense matter equation of state 
(EOS), as the measurement of unprecedentedly high mass 
neutron stars provides support for a subset of "stiff" EOS. 
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Table 4 

Shapiro-delay-derived Mass Measurements 
Source №12 m, (Mo) №15 т, (Мо) c Change? 
J0613—0200 4,5233 3.2114 —0.50 
J0740--6620^ 2.57103 1.99+007 EN 
11600—3053 3,2244 1.207034 -12c 
J1614—2230 1.92210015 1.937100 +1.00 
J1630+3734 ЛЕ 6.3132 e 
J1640+2224 5.2150 2:823 —0.60 
J17414-1351 15297 0.831026 —1.0c 
J1811—2405 2.105. 
11853--1303 eL 1.8478 ... 
В1855+09 1.531+0.098 1.563008 +0.30 
J1903+0327 0.787634 0.944041 +0.50 
J1909—3744 1,5240022 1.574099 +2.20 
J1918—0642 1.241010 1.313009 +0.70 
J1946+3417 3.9104 3.9794 0с 
J2017+0603 22725 2.0773 0c 
J20434-1711 1.641014 1.621010 0.20 
J2302+4442 5.6130 5.2132 0.10 
Notes. 


è №015 т = NG12 m, divided by the average of the NG12 upper and lower 
1c uncertainties. 
P See Section 5.1.4. 


Timed by NANOGrav since 2014, this source is one of the 
most massive neutron stars known. Cromartie et al. (2020) 
reported its mass to be 2.147019 Мо after a series of follow-up 
campaigns around superior conjunction with the GBT. FCP21 
provided an updated measurement of its mass, 2.08*007 Mo, 
after combining the aforementioned follow-up observations 
with a preliminary МС15 and daily cadence Canadian 
Hydrogen Intensity Mapping Experiment pulsar observing 
system (CHIME/Pulsar; CHIME/Pulsar Collaboration et al. 
2021) observations. Radio-derived neutron star mass con- 
straints can be employed as priors in mass-to-radius measure- 
ments from X-ray light-curve modeling (e.g., Riley et al. 2021; 
Miller et al. 2021). 

This work demonstrates that, based on NANOGrav observa- 
tions alone, the mass of PSR J07404-6620 is constrained >7 
times better by NG15 compared to NG12, an improvement that 
can be attributed to the orbital-phase-targeted campaigns first 
presented in Cromartie et al. (2020) and the additional ~3 yr of 
regular NANOGrav observations. However, this analysis does 
not incorporate the high-cadence CHIME data, nor does it take 
into account the newly measured B, (which is dependent on m,, 
Mp, and the distance d to the source) or optimized DM 
modeling as FCP21 does. We therefore do not regard the 
present measurement, which is consistent with the FCP21 to 
within lo, as superseding that result. 

PSR J1614—2230, the first directly observed 2 Mo neutron 
star (Demorest et al. 2010), is timed as part of the NANOGrav 
PTA. In this data release, the addition of newly measured RN 
to its timing model coincides with a ~1o increase in mp. 

Conducting orbital-phase-specific observing campaigns 
around superior conjunction results in demonstrated improve- 
ments in m, constraints for sufficiently highly inclined systems. 
For this reason, five MSPs currently timed by NANOGrav that 
show borderline m, constraints (including PSR J1630+3734, 
PSR J1811—2405, and PSR J1853+ 1303) were subject to such 
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a targeted GBT campaign in spring 2022, the results of which 
will be included in a future manuscript. 


5.1.5. Changes in Noise Parameters since the 12.5 yr Data Set 


Section 43 and the detector characterization paper 
(NG15detchar) describe the WN and RN modeling conducted 
for each pulsar in the data set. A small fraction of WN 
parameters (28/159, 17/159, and 5/159 for EFAC, EQUAD, 
and ECORR, respectively) changed by >30 between NG12 
and NGI5. Even so, these changes are not unexpected as the 
length of the data set increases, especially for pulsars that were 
newly added to NG12. 

Twenty-three pulsars in МО15 were found to have 
significant levels of RN (see filled points in Figure 4). Ten of 
these measurements are newly significant in NG15, while 13 of 
the 14 sources with significant RN in NGI2 continued to favor 
inclusion of RN in their timing models. Only PSR J2317 
+1439, which favored RN amplitude an order of magnitude 
lower than the next-flattest spectral index in NG12, as well as a 
steep spectral index comparable to that of PSR J00304-0451, 
had RN parameters removed from its timing model in NGI5. 
RN parameter values for pulsars with RN in both №012 and 
NGI5 were checked for consistency. No spectral index or 
amplitude values differed between the data sets by more than a 
factor of just over ~3а. Not only are these inconsistencies 
small, but they are also expected to accompany a growing 
data set. 


5.1.6. Kopeikin Parameter Analysis for PSR J1713+0747 


Unlike the majority of the binary systems analyzed here, for 
PSR J1713+0747 we need to use the DDK model that 
incorporates the effects of proper motion and parallax on the 
binary orbit (Kopeikin 1995, 1996; see also Appendix D). 
However, in addition to the convention ambiguity discussed in 
Appendix D, the nature of the analysis in Kopeikin (1996) also 
leads to further ambiguity with nearly similar solutions possible 
for different choices of inclination i and longitude of the 
ascending node €). These are locations where the change in 
projected semimajor axis óx is equal to the best-fit value, but 
where 6x is different. The locations of these points can be found 
by first determining where óx is 0, which is where 


(5) 


(defined in ecliptic coordinates). With an initial best-fit pair (i, 
О) the other points can then be identified as 


i Q 
i 2% – 9 — 180° 
— 
" 1800— i — Q4 180? (6) 
i80 —i 20-0 


We wished to verify that the solution we had settled on was in 
fact a global and not just a local minimum. To do this, we 
computed a grid of ж over values of (i, О) with PINT, with all 
other parameters freely fit at each grid point. We used a 
preliminary version of the narrowband TOAs for this purpose, 
but subsequent changes and comparison with wideband TOAs 
did not indicate any difference. 
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Figure 4. Comparisons of the two RN parameters (the spectral index 7,.4 and 10210 amplitude, 102, Area) between NG15 and NGI2. Гей: logi) Arq VS- Yrea for NG15 
(blue) and NG12 (orange) without error bars. Filled squares are values for pulsars whose RN was found to be significant and is included in the timing model. Crosses 
are the measured RN parameters for pulsars without significant RN. Middle and right panels: same as the left panel, but split into the NG15 values (middle) and NG12 
values (right), including the 10 confidence interval error bars. For crosses (insignificant RN) appearing in the “high-A,.q” group of points, the RN priors are the only 
information constraining those measurements, that is, the data provide no information about the RN content of those MSPs. 


We stepped through all values of 2 from —180° to +180° in 
1° increments and similarly stepped through values of i from 0° 
to 180° in 1° increments. As mentioned in Appendix D, our 
final fit with PINT used ecliptic coordinates, but we wished to 
compare to results from the literature that were done in 
equatorial coordinates. The value of i will not change, but the 
value of €) will change. However, at this location in the sky the 
equatorial (ICRS) and ecliptic coordinate systems are almost 
aligned, with OQicrs = Осирис = 523. 

The results of this fit are shown in Figure 5. Тһе four local 
minima given by Equation (6) are clearly visible, and we 
actually found that our initial fit had ended up in the wrong 
local minimum. However, we have now verified that the 
minimum identified by our nonlinear fitter is the global 
minimum, with x? differences of 261-904 for the alternatives 
(with approximately 59,000 degrees of freedom). 

This also agrees with the solutions from Splaver et al. (2005) 
and Alam et al. (2021a), once they have been corrected to use 
the Damour & Taylor (1992, hereafter DT92) convention for 
defining О (Appendix D) and ecliptic coordinates, again 
verifying that we have found the global minimum. We note that 
the Alam et al. (2021a) result was derived using the T2 timing 
model and that the resulting values of О and i may be shifted 
from the values found by the other models by roughly 0.12, 
due to the use of unadjusted Keplerian parameters in 
computing binary delays (see Appendix D). Finally, we also 
get the same result if we fit using Tempo and convert from the 
IAU to DT92 convention for defining Q (Appendix D). 


5.2. Comparison of Split-telescope Data Sets 


Most pulsars in the data set are observed by only one 
observatory, either GBT or Arecibo. However, two pulsars are 
observed by both GBT and Arecibo (PSR B19374-21 and PSR 
J1713--0747), and several are observed with the VLA in 
addition to GBT or Arecibo. 

For pulsars observed with GBT and Arecibo, we create 
timing solutions using TOAs from each telescope separately, in 
addition to the combined data set, and compare respective 
parameters. These data sets enable separate GBT and Arecibo 
"split-telescope" GWB searches, as in NGl5gwb. Equally 
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Figure 5. Grid of x? as a function of inclination i and longitude of the 
ascending node О for PSR J17134-0747, computed with PINT. АП other 
parameters were fit at each grid point. The values are in the DT92 (Damour & 
Taylor 1992) convention (see Appendix D). The color scale and contours are 
the logarithm of difference of X^ compared to the best-fit value, with an inset 
zooming in to the best-fit region. The magenta diamond is the value from the 
12.5 yr data release (Alam et al. 2021a) converted from the IAU convention 
to DT92; these parameters will be slightly shifted owing to the choice of timing 
model. The cyan pentagon shows the value from Splaver et al. (2005) 
converted from the IAU convention to DT92 and additionally converted from 
equatorial to ecliptic coordinates. The green hexagon is the best fit when using 
equatorial coordinates, converted to ecliptic. The black square is the best fit 
when using ecliptic coordinates directly. Finally, the blue star is the minimum 
point in the x? grid. 


important, we use these separate data sets to check consistency 
between observatories. As both telescopes have extensive data 
sets for PSR B1937--21 and PSR J1713--0747, we expect the 
solutions to be statistically consistent with each other and with 
the combined solution. Comparisons of these solutions have 
confirmed that this is the case. 
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Figure 6. Geometry of the binary system, showing the position angle of ascending node in the ТАЏ and DT92 conventions, Оду and Орто2; the inclination angle in 
the ТАЏ and DT92 conventions, irau and 1792; and the angle of periastron, w, which is the same in both conventions. The cardinal directions—N, S, E, and W for 
north, south, east, and west— are relative to whichever astrometric coordinate system is used for the pulsar position fit (equatorial or ecliptic). Red portions of the 
figure are in the plane of the sky; blue portions are in the orbital plane. This figure is adapted from Splaver et al. (2005). 


Additionally, for pulsars observed with the VLA and 
Arecibo (J1903+0327) or the VLA and GBT (Ј1600–3053, 
J1643—1224, and J1909—3744), we create separate Arecibo or 
GBT-only data sets. Due to the small number of TOAs, we do 
not create VLA-only data sets for these sources. We again 
compare the Arecibo or GBT-only solution to the complete 
solution. In all cases, we find that parameters are consistent 
within 3c, except where differences in the systems used, data 
spans, frequency ranges, or RN parameters would make 
changes expected. 


6. New Pulsars 


Since the 12.5 yr data set, many new MSPs are now being 
observed, and results for 21 new pulsars are included in NG15. 
Historically, NANOGrav has had a close relationship with 
pulsar survey groups using the Arecibo and Green Bank 
Observatories, and іп NGI5 new sources came from three 
untargeted surveys at these telescopes, the Arecibo 327 MHz 
Drift-Scan Pulsar Survey (PSR J05094-0856 and PSR J0709 
+0458; Martinez et al. 2019), the Pulsar Arecibo L-band Feed 
Array survey (PSR J05574-1551; Scholz et al. 2015), and the 
Green Bank North Celestial Cap survey (PSR J04064-3039). 
Four more came from targeted surveys at Arecibo/Green Bank, 
guided by Fermi Large Area Telescope (LAT) unassociated ^- 
ray sources (PSR J06054-3757, PSR J13124-0051, PSR J1630 
+3734, and PSR J0614—3329; Ransom et al. 2011; Sanpa- 
Arsa 2016). We also draw on discoveries from similar targeted 
surveys carried out with the Parkes Telescope (PSR Ј1012 
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—4235; Camilo et al. 2015) and the Effelsberg Telescope 
(PSR J1745+1017; Barr et al. 2013) and additional blind 
surveys like the High Time Resolution Universe Survey (PSR 
J1705—1903, PSR J1719—1438, and PSR J1811—2405; Ng 
et al. 2014; Morello et al. 2019) and the Parkes Multibeam 
Pulsar Survey (PSR J1802—2124; Faulkner et al. 2004), which 
both use the Parkes Telescope. Some sources added here are 
also monitored by the EPTA (PSR J1751—2857 and PSR 
J1843—1113; Desvignes et al. 2016) and the PPTA (PSR J0437 
—4715; Reardon et al. 2021). PSR J0610—2100, PSR 11022 
+1001, PSR J1730—2304, and PSR J2124—3358 are also 
monitored by both EPTA and PPTA. While seemingly 
redundant, overlap in observing programs among regional 
PTAs like this is useful for International PTA (IPTA) data 
combination and for diagnosing related issues (see, e.g., Hobbs 
et al. 2010; Perera et al. 2019). Overlap in observing programs 
can enhance our frequency coverage and cadence for these 
pulsars as well. 


7. Summary and Conclusions 


In this paper we have presented the 15 yr data set from the 
NANOGrav Collaboration, using data from the Green Bank, 
Arecibo, and VLA telescopes. Our longest timing baselines 
have increased by ~3 yr over our previous data set. Even more 
impactful is the increase in the number of pulsars, from 47 
MSPs in the 12.5 yr data set to 68 MSPs in the present data set; 
this is the largest fractional increase (~50%) in the number of 
pulsars since the size of the array doubled between NG5 and 
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NG9. This large increase was made possible by NANOGrav’s 
push to evaluate many new MSPs discovered in large-area sky 
surveys and targeted gamma-ray-guided surveys; it is notable 
because a PTA’s sensitivity to a GWB increases linearly with 
the number of pulsars in the array and as the square root of the 
overall timing baseline (Siemens et al. 2013). 

Another high-impact addition in this data set is the novel 
pipeline, PINT Pal, that automates the timing process. The 
user interfaces with a configuration file to set up the run by 
specifying which TOA files, noise chains, etc., to use, and the 
pipeline runs the timing analysis and outputs a series of 
diagnostic plots and tables with which the user determines 
whether or not the full timing model is sufficient. The timing 
analysis can easily be rerun by editing the configuration file, 
e.g., to conduct iterative noise analyses when parameters are 
added or removed from a given pulsar’s timing model. This 
pipeline lowers the bar of entry for students and others who are 
relatively new to pulsar timing, and it will be of great use in 
future data sets, especially as the number of pulsars continues 
to grow. The intent is for this pipeline to also be adapted and 
used by other PTA collaborations and for more general pulsar 
timing analysis and IPTA data combination. The publicly 
available PINT_Pal can be found at https://github.com/ 
nanograv /pint_pal. 

In NG12, we included a comparison between the Tempo and 
PINT timing results. A thorough comparison between timing 
software was also included in Luo et al. (2019). In the 15 yr 
data set, we have shifted to using PINT as the primary timing 
software. PINT is a modular, Python-based software, and this 
transition decreases barriers to entry in understanding our 
pipeline and results. In NG12, we conducted systematic 
comparisons between Tempo and PINT, demonstrating the 
consistency of our prior and current timing software. 

In Section 5, we performed а number of comparisons: we 
compared the astrometric, binary, and RN parameters measured 
in the 12.5 and 15 yr data sets and the timing and RN parameter 
results between individual-telescope data sets for those pulsars 
observed with both Arecibo and GBT. Overall, the parameters 
are consistent between the various data sets. A small number of 
parallax, proper-motion, and mass measurements have changed 
significantly between the 12.5 and 15 yr data releases, as 
detailed in Sections 5.1.1 and 5.1.3. In МСІ2, we reported 
detectable levels of RN for 14 pulsars; in the 15 yr data set we 
measure RN for 23 pulsars, 13 of which were also found to 
show RN in NG12. Aside from the one source for which we 
detected RN in NG12 but not in the 15 yr data set, RN 
parameters were consistent between the data sets. 

NANOGrav continues to be committed to public data 
releases for GW detection, as well as for individual pulsar 
studies." In addition, we are committed to data sharing and 
collaborative analysis within the IPTA, of which the EPTA, 
InPTA, and PPTA are also members. NGI5, along with the 
most recent EPTA, InPTA, and PPTA data sets, have already 
been shared with the IPTA in order to create a combined data 
set, which will be more sensitive to GW signals than any 
individual РТА data set. IPTA data combination analyses have 
begun. 

We are also committed to continuing to increase the size of 
the NANOGrav PTA. The loss of Arecibo was significant, as 
new faint pulsars that fall in its decl. range can no longer be 


55 Data from this paper are available at ћир: / /data.nanograv.org and preserved 
on Zenodo at doi:10.5281/zenodo.7967585. 
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added to the array, and we have been unable to continue 
observing a small number of faint pulsars previously observed 
at Arecibo. Despite this, NANOGrav will continue to evaluate 
new MSPs for inclusion in our GBT, VLA, and CHIME timing 
campaigns. With the proposed DSA-2000 telescope (Hallinan 
et al. 2019), NANOGrav will further expand its timing 
campaign, increasing the number of MSPs observed and thus 
our sensitivity to GWs. 
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noise modeling results. N.S.P. undertook comparisons 
between PINT, Tempo, and Tempo2 results. D.L.K. and I. 
H.S. aided in comparisons specifically related to Kopeikin 
parameters, as described in Section 5.1.6. M.E.D. performed 
the comparison of astrometric parameters between data sets. 
H.T.C., K.C., and J.K.S. performed comparisons between the 
narrowband and wideband data sets. H.T.C. and E.F. 
undertook analysis of binary systems. D.C.G. compared 
split-telescope results for the subset of pulsars observed with 
multiple telescopes. Significant contributors to the text, 
tables, and figures in this paper include H.T.C., M.E.D., P. 
B.D., D.C.G., J.G., J.S.H., D.L.K., M. T.L., M.A.M., D.J.N., 
S.M.R., LH.S., and J.K.S. H.T.C. produced the timing 
residual plots. J.G. and N.G. oversaw and maintained much 
of the computational infrastructure related and essential to 
this work and contributed to the design and public release of 
the PINT_PAL code repository. We dedicate this paper to our 
wonderful colleague Jing Luo, whom we lost this past year. 
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Appendix A 
Further Details on TOA Removal 


А.1. Poor Front-end/Back-end Cut (poorfebe) 


In most cases, the poorfebe cut was applied in cases 
where the vast majority of, or all, TOAs from a given front- 
end/back-end combination were being cut by other means. It 
was only used for six MSPs and mostly to indicate poor 
performance at higher observing frequencies (e.g., 3/6 GHz 
with the VLA and 2 GHz with Arecibo). In one case, a pulsar 
had been observed more than three times with a nonstandard 
front end owing to RFI issues at Arecibo Observatory and 
ad hoc scheduling changes. These TOAs were not caught by 
our orphan data cut, and for consistency we removed them with 
poorfebe. 


A.2. Bad Range Cut (badrange) 


In preliminary stages of our analysis, we discovered that a 
significant portion of data collected at Arecibo exclusively 
between 2017 August and 2018 November were corrupted, and 
the issue was traced back to a malfunctioning local oscillator 
(LO).*^ During this time period, (һе LO's reference frequency 
exhibited sudden shifts by 5-10 MHz, and sometimes erratic 
wandering across this range on submillisecond timescales. This 
sort of behavior manifests itself as what looks like improperly 
dedispersed archives in the NANOGrav data set and causes 
often unpredictable amounts of smearing and phase variations 
within individual observations. To determine the scope of the 
problem, we inspected phase versus frequency residuals, after 
subtracting the pulse portrait from the corresponding band. 
Corrupted scans were easily detectable by eye for our brightest 
pulsars (i.e., PSR J17134-0747), but due to some concern about 
low-level effects on high-precision timing for all of our pulsars 
monitored at Arecibo, we took a conservative approach and 
excised all scans within the time range MJD 57984—58447, 
where corruption was apparent. TOAs from this MJD range 
were assigned the cut flag badrange. After identifying this 
issue in МС15, we have implemented additional quality 
assurance measures to ensure that we catch similar issues 
more quickly in the future. 


A.3. Orbital Phase Range Cut (eclipsing) 


Several black widow pulsars included in NANOGrav's 
regular observing schedule appear in NG15, and two of these, 
PSR J1705—1903 and PSR J1802—2124, exhibit eclipses. By 
examining timing residuals as a function of orbital phase and 
looking at individual scans, we imposed an eclipsing cut so 
that TOAs generated within 10956-1596 of an orbit from 
superior conjunction (when the pulsar's signal is eclipsed by its 
companion) are automatically removed. 


А.4. Gibbs Outlier Cut (outlier10, maxout, epochdrop) 


While NG12 employed the Vallisneri & van Haasteren 
(2017) outlier identification algorithm, the growing TOA 
volume of our narrowband data set, as well as the desire to 


59 Diagnosing this issue was aided by private communication with Shriharsh 
Tendulkar and others. 
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iteratively perform outlier analyses in response to timing model 
modifications, necessitated a more computationally efficient 
approach. In this work, we use a Gibbs sampler to determine 
TOA outliers in a fully Bayesian manner that is more 
computationally efficient than acceptance /rejection-based sam- 
plers. Wang & Taylor (2022) provide a thorough overview of 
this technique, including a demonstration of its efficacy when 
applied to the NANOGrav data set. Both outlier methods 
(Vallisneri & van Haasteren 2017; Wang & Taylor 2022) are 
available in enterprise outliers," which is a depen- 
dency for our timing analysis pipeline (see Section 4 for more 
details). As in NG12, any narrowband TOAs with outlier 
probability p > 0.1 were removed and assigned the cut flag 
outlier10. Through experimentation, we found that in situ 
ations where more than five narrowband TOAs were flagged as 
outliers, all TOAs from that file were often corrupted. 
Therefore, if more than 8% of TOAs (usually 5/64, but there 
are fewer total TOAs generated for some bands) from a single 
observing epoch were flagged for removal by the Gibbs 
algorithm, the remaining TOAs were also cut using the 
maxout flag. Neither of these steps was used for our 
wideband data set, but since there are one to two wideband 
TOAs generated per file, the epoch F-test cut (see NGI2, 
Section 2.5.8) provides a similar per-TOA outlier assessment in 
that case (epochdrop). This final stage of our automated 
outlier analysis is also applied to narrowband TOAs. 


А.5. Bad TOA/File Cut (badtoa, badfile) 


After the aforementioned outlier removal techniques were 
applied and the timing model was fully fit, a comparatively 
small number of TOAs/files were removed manually with 
badtoa and badfile cut flags, respectively. For every 
TOA/file removed this way, corresponding data cubes were 
inspected and the exact TOA/file and reason for removal were 
recorded in timing configuration files for transparency and 
posterity. Reasons for manual TOA/file cuts include the 
following: snr/non-detection when TOAs barely 
exceeded the S/N threshold and/or where no signal was 
visible by eye; rfi when RFI was still clearly present; 
few chans when only one or several channels remained 
post-zapping, suggesting that almost the entire band had been 
affected by RFI; smeared when archives had been folded 
improperly or significant pulse broadening was apparent for 
other reasons; and isolated when one or several observa- 
tions were separated from the rest by a long time span (e.g., test 
observations for PSR J1730—2304 were separated from the 
regular monitoring by 12 yr). A total of 14 narrowband and 22 
wideband TOAs were removed despite there being no obvious 
reason for removal (unknown), but these have also been 
explicitly recorded in the configuration files. 


Appendix B 


This appendix includes summaries of basic pulsar para- 
meters and narrowband TOA statistics (Table 5) and both 
narrowband and widebandtiming model fits (Table 6). 


37 https: / [ github.com/nanograv/enterprise outliers 
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Source 


10023--0923 
10030--0451 
10340--4130 
10406--3039 
10437-4715 
10509--0856 
10557--1551 
10605--3757 
10610-2100 
10613-0200 
J0614—3329 
0636-5128 
0645-5158 
J0709--0458 
10740--6620 
10931-1902 
J1012+5307 
J1012—4235 
J1022+1001 
J1024—0719 
J1125+7819 
J1312+0051 
J1453+1902 
J1455—3330 
J1600—3053 
J1614—2230 
J1630+3734 
J1640+2224 
J1643—1224 
J1705—1903 
J1713+0747 
J1719—1438 
J1730—2304 
11738--0333 
J17414-1351 
J1744—1134 
J1745+1017 
J1747—4036 
11751-2857 
1802—2124 
1811—2405 
11832—0836 
11843—1113 
J1853+1303 
B1855+09 

11903--0327 
J1909—3744 
J1910+1256 
J1911+1347 
J1918—0642 
11923--2515 
В1937--21 

11944--0907 
11946--3417 
В1953--29 

Ј2010— 1323 
J20174-0603 
J2033+1734 
J2043+1711 
J2124—3358 


P 
(ms) 


3.05 
4.87 
3.30 
2.61 
5.76 
4.06 
2.56 
2.73 
3.86 
3.06 
3.15 
2.87 
8.85 
34.43 
2.89 
4.64 
5.26 
3.10 
16.45 
5.16 
4.20 
4.23 
5.79 
7.99 
3.60 
3.15 
3.32 
3.16 
4.62 
2.48 
4.57 
5.79 
8.12 
5.85 
3.75 
4.07 
2.65 
1.65 
3.91 
12.65 
2.66 
2,12 
1.85 
4.09 
5.36 
2.15 
2.95 
4.98 
4.63 
7.65 
3.79 
1.56 
5.19 
3.17 
6.13 
5.22 
2.90 
5.95 
2.38 
4.93 


dP/dt 
(10729) 


1.14 
1.02 
0.70 
0.83 
5.73 
0.44 
0.72 
0.47 
1.23 
0.96 
1.74 
0.34 
0.49 
38.02 
1.22 
0.36 
1.71 
0.66 
4.33 
1.86 
0.69 
1.75 
1.17 
2.43 
0.95 
0.96 
1.07 
0.28 
1.85 
2.15 
0.85 
0.80 
2.02 
2.41 
3.02 
0.89 
0.26 
1.31 
1.12 
7.29 
1.34 
0.83 
0.96 
0.87 
1.78 
1.88 
1.40 
0.97 
1.69 
2:57 
0.96 
10.51 
1.73 
0.31 
2.97 
0.48 
0.80 
112 
0.52 
2.06 


DM 
(pc ст?) 


14.4 
4.4 
49.6 
49.4 
2.6 
38.3 
102.6 
20.9 
61.3 
38.8 
37.1 
11.1 
18.2 
44.3 
15.0 
41.5 
8.9 
717 
9.4 
8.4 
11.2 
15.3 
14.1 
13.6 
52.3 
34.5 
14.1 
18.5 
62.3 
57.5 
16.0 
36.8 
9.6 
33.8 
24.2 
3.1 
24.0 
153.0 
42.8 
149.6 
60.6 
28.2 
60.0 
30.6 
13.3 
297.6 
10.4 
38.1 
31.0 
26.6 
18.9 
71.1 
24.4 
110.2 
104.5 
22.2 
23.9 
25.1 
20.7 
4.6 


РЬ 
(days) 


0.1 


7.0 
5.7 
4.9 
4.8 
55.7 
0.3 
1.2 
53.6 
0.1 
4.4 
4.8 
0.6 
38.0 
7.8 
15.4 
38.5 
76.2 
14.3 
8.7 
12.5 
175.5 
147.0 
0.2 
67.8 
0.1 
0.4 
16.3 
0.7 
110.7 
0.7 
6.3 


115.7 
12.3 
95.2 

1.5 

58.5 


10.9 
27.0 
117.3 
2.2 
56.3 
1.5 


327 MHz 


Table 5 
Basic Pulsar Parameters and Narrowband TOA Statistics 


Median Scaled TOA Uncertainty? (us)/ Number of Epochs 


430 MHz 


0.067/81 
0.188/257 


0.140/15 


1.211/54 


0.057/264 


0.154/90 


0.369/90 
0.227/136 


0.635/60 
0.282/73 
0.307 /80 
0.321/78 
0.218/6 
0.233/59 
0.096/205 
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820 MHz 


0.888/102 


1.205/20 
0.894/37 
0.111/165 
0.554/29 
0.269/72 
0.299/98 
0.489 /127 
1.020/84 


0.406/169 
1.129/18 


0.644/120 
1.118/75 


1.224/149 
0.285/144 
0.397/127 
0.586/29 


0.306/161 
0.288/32 
0.195/161 
0.958/35 
0.561/37 


0.161/158 
1.092/94 
1.825/25 
0.955/35 
0.375/38 


0.605 /85 
0.578/38 


0.072/154 
0.533/155 
0.008/161 
0.430/125 


0.809 /35 


1.4 GHz 


0.585/91 
0.491 /285 
2.240/101 
0.416/38 
0.081/20 
1.319/33 
1.402/40 
1.778/23 
1.319/36 
0.654/165 
1.044/26 
0.596/73 
0.823/100 
3.105/50 
0.792/165 
1.840/83 
0.781/175 
2.081/27 
0.428/50 
1.132/124 
275/73 
2.012/45 
2.426/68 
2.100/145 
0.253/148 
0.677/142 
1.080/33 
0.428/283 
0.567/164 
0.192/31 
0.089/676 
1474/31 
1.187/36 
0.563/99 
0.338/112 
0.258/159 
0.703/47 
1.285/93 
1.824/38 
0.978/33 
0.657/37 
0.585/88 
0.580/39 
0.691 /97 
0.241/146 
0.561/98 
0.131/390 
0.406/105 
0.192/68 
1.000/162 
1.055/92 
0.015/279 
0.991/103 
0.507/67 
0.948/97 
1.069 /123 
0.464/85 
1.377/67 
0.299/227 
1.656/40 


2.1 GHz 


1.171/94 
0.918 /30 


4.766/31 
1464/14 


7.486/40 


0.605/35 


2.319/32 
7.343/1 


0.080/263 
0.776/77 
0.139/10 


1.498 /43 


0.591 /96 


0.791/96 


0.022/98 
1.568/12 
0.730/60 
2.076/5 


0.595 /60 
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3.0 GHz 


0.080/28 


0.966 /23 


2.012/23 


0.364/27 


2.004/12 
0.226/31 


Span 
(yr) 


9.0 
15.5 
8.1 
3.6 
4.8 
3.6 
4.6 
34 
34 
15.0 
24 
6.3 
8.9 
4.6 
6.3 
7.1 
15.5 
34 
5.6 
10.5 
6.3 
4.6 
7.0 
15.7 
12.5 
11.5 
3.5 
15.5 
15.7 
3.7 
15.5 
34 
34 
10.7 
11.0 
15.7 
4.5 
8.1 
3.5 
3.5 
3.5 
7.1 
3.5 
9.1 
15.6 
10.7 
15.5 
114 
7.0 
15.5 
9.0 
15.9 
12.5 
57 
11.1 
10.5 
8.3 
7.0 
9.1 
3.5 
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Table 5 
(Continued) 

Source Р dP/dt DM P, Median Scaled TOA Uncertainty? (из) /МатБег of Epochs Span 

(ms) (10729) _ (рест ?) (days) 327 MHz 430 MHz 820 MHz 1.4 GHz 2.1 GHz 3.0 GHz (yr) 
J2145—0750 16.05 2.98 8.9 6.8 0.322/145 0.688/145 tee 15.5 
J2214+3000 3.12 1.47 22.5 0.4 Ut tee 0.844/91 1.203/61 8.8 
12229--2643 2.98 0.15 22.7 93.0 0.280/69 1.273/76 Ut 7.0 
12234--0611 3.58 1.20 10.8 32.0 0.449/55 0.271/65 ves 6.5 
12234--0944 3.63 2.01 17.8 0.4 © Ut 0.444/62 0.840/60 7.1 
J2302+4442 5.19 1.39 13.8 125.9 vee 1.200/101 2.569/97 © 8.1 
J2317+1439 3.45 0.24 21.9 2.5 0.084/71 0.079/274 0.714/241 Ut 15.6 
J2322--2057 4.81 0.97 13.4 vee 0.356/56 1.148/58 1.481/8 5.4 

Nominal scaling factors” for ASP/GASP: 0.58 0.45 0.80 0.80 0.80 
GUPPI/PUPPI: 0.71 0.49 1.34 2.49 2.14 КЕ 
YUPPI: E 2.83 4.12 

Notes. 


а Original narrowband TOA uncertainties were scaled by their bandwidth-time product | 


bandwidths and integration times. 
> TOA uncertainties can be rescaled to the nominal full instrumental bandwidth by dividing by these scaling factors. 


Av z 
100 MHz !800 


1/2 
) to remove variation due to different instrument 


Table 6 
Summary of Timing Model Fits* 
Somme Number Number of Fit Parameters” rms° (us) Red Noise? Figure 
of TOAs 5 А В DM FD J Full White Аса "Уеа logioB Number 
10023--0923 15,896 3 5 10 92 4 1 0.326 0.24 7 
824 3 5 10 89 0 1/2 0.320 —0.11 
10030--0451 19,579 3 5 0 289 4 2 0.856 0.251 0.003 —4.7 >2 8 
727 3 5 0 289 0 2/3 0.794 0.263 0.003 —4.9 >2 
10340+4130 11,093 3 5 0 108 2 1 0.597 —0.22 9 
228 3 5 0 108 0 1/2 0.591 —0.12 
J0406+3039 2446 3 5 5 39 1 1 0.176 —0.09 10 
71 3 5 5 39 0 1/2 0.300 —0.21 
J0437—4715 5830 3 5 T 30 5 1 0.186 0.110 0.027 —0.1 >2 11 
117 3 5 6 31 1 1/2 0.123 0.27 
10509--0856 2169 3 5 5 38 0 1 0.695 0.73 12 
66 3 5 5 37 0 1/2 0.760 0.17 
J0557--1551 525 3 5 5 41 0 1 0.353 0.01 13 
47 3 5 5 43 0 1/2 0.233 —0.13 
10605--3757 554 3 5 5 26 0 1 1.153 —0.11 14 
47 3 5 5 29 0 1/2 1.125 —0.08 
J0610—2100 4885 3 5 5 38 4 1 1.703 1.059 0.229 —2.6 >2 15 
214 3 5 5 38 2 1/2 2.162 1.025 0.105 —4.5 >2 
10613—0200 17,124 3 5 8 174 2 1 0.749 0.168 0.022 —2.9 22 16 
423 3 5 8 174 0 1/2 0.704 0.151 0.013 —3.2 >2 
0614—3329 1714 3 5 5 30 1 1 0.276 —0.08 17 
56 3 5 5 30 0 1/2 0.350 —0.06 
J0636+5128 32,222 3 5 8 77 1 1 0.664 —0.14 18 
1221 3 5 8 79 0 1/2 0.645 0.24 
J0645+5158 17,670 3 5 0 113 2 1 1.592 0.49 19 
289 3 5 0 114 0 1/2 0.164 1.55 
10709--0458 3030 3 5 5 51 1 1 1.165 0.01 20 
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Table 6 
(Continued) 
Sonic Number Number of Fit Parameters” rms^ (us) Red Noise? Figure 
of TOAs S A B DM FD J Full White Аса “гей logioB Number 
91 3 5 5 52 0 1/2 1.151 et e e 0.00 
10740--6620 13,401 3 5 8 109 3 1 0.286 e E e —0.18 21 
360 3 5 8 109 0 1/2 0.296 et E e —0.22 
J0931—1902 5473 3 5 0 91 0 1 0.415 et E e. —0.29 22 
190 3 5 0 91 0 1/2 0.382 e Ue e —0.11 
J10124-5307 25,837 3 5 7 177 4 1 0.925 0.289 0.219 —0.9 22 23 
628 3 5 7 171 3 1/2 0.967 0.224 0.241 —0.8 >2 
J1012—4235 797 3 5 5 28 0 1 0.708 ... E et —0.09 24 
65 3 5 7 35 0 1/2 0.623 e E e —0.11 
J1022+1001 3978 3 5 6 55 5 2 2.835 e E e 0.15 23 
116 3 5 6 56 0 2/3 3.050 ЕЕ E e. 0.09 
J1024—0719 12,635 3 5 1 134 5 1 0.239 te E et —0.16 26 
288 3 5 1 134 0 1/2 0.245 e. E e. —0.22 
J11254-7819 8723 3 5 6 78 2 1 0.688 te e e 0.10 27 
206 3 5 6 79 0 1/2 0.619 e E e. —0.11 
J13124-0051 1705 3 5 5 48 0 1 0.731 et E e. —0.14 28 
76 3 5 5 49 0 1/2 0.632 e. -- e —0.11 
J14534-1902 2551 3 5 0 69 0 1 0.783 et E e. —0.18 29 
122 3 5 0 71 0 1/2 0.895 e -- er —0.09 
J1455—3330 10,818 3 5 6 157 1 1 0.735 e E e —0.14 30 
357 3 5 6 156 0 1/2 0.663 et E e —0.17 
J1600—3053 22,955 3 5 9 181 3 2 0.271 0.202 0.037 —0.9 22 31 
481 3 5 8 182 1 2/3 0.338 0.145 0.052 —1.3 >2 
1614—2230 18,445 3 5 8 150 0 1 0.354 0.211 0.010 —3.0 >2 22 
367 3 5 8 151 0 1/2 0.311 0.202 0.010 -2.6 >2 
J1630+3734 1815 3 5 8 34 1 1 0.271 et E e. —0.09 33 
71 3 5 7 38 0 1/2 0.473 et E e. —0.04 
J1640+2224 14,066 3 5 9 284 4 1 0.200 e E e —0.06 34 
609 3 5 8 283 0 1/2 0.181 e. E e —0.16 
J1643—1224 22,144 3 5 7 193 5 2 2:335 0.898 0.543 —0.7 >2 35 
478 3 5 6 186 1 2/3 2.598 0.675 0.538 —1.2 >2 
11705—1903 9871 3 5 10 43 2 1 1.124 0.226 0.217 —0.4 >2 36 
253 3 5 10 45 0 1/2 1.231 0.252 0.296 —0.8 22 
J17134-0747 59,389 3 5 8 398 4 5 0.201 0.095 0.011 —2.2 >2 37 
1495 3 5 8 398 3 5/6 0.199 0.093 0.009 —2.8 >2 
1719—1438 6356 3 5 5 42 1 1 2.471 e -- e —0.10 38 
463 3 5 6 41 0 1/2 2.703 e E e. —0.08 
J1730—2304 4870 3 5 0 44 2 1 0.277 e E e 0.61 39 
93 3 5 0 43 1 1/2 0.310 e e e —0.01 
J17384-0333 8790 3 5 5 104 1 1 0.262 te e e 1.57 40 
336 3 5 5 102 0 1/2 1.022 0.259 0.000 —6.9 22 
J17414-1351 5582 3 5 9 113 2 2 0.182 e E e. 0.13 41 
208 3 5 9 104 0 1/2 0.169 te E -- 0.13 
J1744—1134 17,745 3 5 0 169 4 1 0.505 0.287 0.022 —2.5 >2 42 
433 3 5 0 169 0 1/2 1.042 0.271 0.006 —4.7 >2 
J1745+1017 3017 3 5 5 49 1 1 12.519 0.850 1.167 —2.5 22 43 
91 3 5 5 49 0 1/2 9.907 0.379 1.095 -2.4 >2 
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Table 6 
(Continued) 
Soure Number Number of Fit Parameters” rms° (us) Red Noise? Figure 
of TOAs S A B DM FD J Full White Área red logioB Number 

J1747 —4036 11,055 3 5 0 106 1 1 3.214 1.480 0.268 -2.4 >2 44 
222 3 5 0 105 0 1/2 2.802 0.756 0.364 —2.0 >2 

J1751—2857 2025 3 5 5 44 1 1 0.560 —0.11 45 
89 3 5 5 45 0 1/2 0.530 —0.06 

J1802—2124 6796 3 5 5 45 1 1 2.247 0.878 0.544 —1.5 >2 46 
126 3 5 7 45 5 1/2 0.974 e 1.09 

J1811—2405 5266 3 5 Т 46 1 1 0.220 -0.09 47 
103 3 5 7 46 0 1/2 0.223 -0.10 

11832-0836 7739 3 9 0 93 1 1 0.214 0.07 48 
207 3 5 0 93 1 1/2 0.208 -0.07 

1843-1113 4595 3 5 0 40 1 1 0.220 -0.09 49 
103 3 5 0 40 0 1/2 0.235 0.12 

11853--1303 4570 3 5 8 98 0 1 0.337 0.163 0.034 -2.0 >2 50 
184 3 5 8 96 0 1/2 0.354 0.171 0.088 —5.4 >2 

B1855+09 7758 3 5 7 147 3 1 0.829 0.330 0.011 —3.7 22 51 
364 3 5 7 149 0 1/2 0.845 0.352 0.016 —3.4 >2 

11903--0327 6856 3 9 8 114 3 2 3.737 0.711 0.627 —1.4 >2 52 
226 3 5 8 114 0 2/3 1.670 0.207 0.318 —1.1 >2 

J1909—3744 35,037 3 5 9 325 1 2 0.303 0.066 0.002 —44 22 53 
833 3 5 9 327 0 2/3 0.287 0.062 0.003 —42 >2 

J1910+1256 6486 3 5 6 114 2 1 0.413 0.06 54 
216 3 5 6 114 0 1/2 0.438 —0.16 

11911+1347 3786 3 5 0 69 2 1 0.074 —0.19 55 
126 3 5 0 69 0 1/2 0.088 —0.26 

J1918—0642 18,875 3 5 8 166 1 1 0.338 e E .. 1.67 56 
487 3 5 7 168 0 1/2 0.398 0.296 0.032 —2.0 >2 

J1923+2515 4001 3 5 0 90 1 1 0.280 0.07 57 
170 3 5 0 92 0 1/2 0.214 —0.25 

B1937+21 23,023 3 5 0 262 5 3 5.774 5.774 0.026 —42 5:2, 58 
660 3 9 0 267 2 3/4 5.958 0.057 0.037 —3.6 22 

J19444-0907 5328 3 5 0 104 2 2 0.461 0.47 59 
180 3 5 0 95 0 1/2 0.411 —0.03 

J1946+3417 4743 3 5 8 72 1 1 1.387 0.305 0.286 —1.4 22 60 
129 3 5 8 71 0 1/2 0.740 0.162 0.154 —1.8 >2 

B1953+29 5126 3 5 6 98 2 2 1.158 0.378 0.216 —1.4 >2 61 
173 3 5 6 94 0 1/2 1.450 0.279 0.195 —2.2 >2 

12010—1323 17,077 3 5 0 142 1 1 0.274 0.01 62 
350 3 5 0 143 0 1/2 0.271 —0.07 

J2017+0603 3512 3 5 T 92 0 2 0.109 —0.17 63 
154 3 5 7 93 0 2/3 0.126 —0.22 

J2033+1734 3847 3 5 5 68 2 1 0.468 —0.09 64 
133 3 5 3 70 0 1/2 0.399 —0.22 

J2043+1711 7398 3 5 8 228 2 1 0.115 0.59 65 
459 3 5 8 223 0 1/2 0.103 0.75 

J2124—3358 4982 3 5 0 40 1 1 0.338 —0.12 66 
104 3 5 0 39 0 1/2 0.531 —0.17 

J2145—0750 18,675 3 5 Т 161 4 1 0.799 0.644 0.111 —0.5 >2 67 
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Table 6 
(Continued) 
Source Nuniber Number of Fit Parameters? rms° (us) Red Noise? Figure 
of TOAs S A B DM FD J Full White Área red logioB Number 

400 3 5 Т. 161 0 1/2 1.045 0.358 0.068 —2.4 >2 

J2214+3000 7425 3 5 5 96 2 1 0.407 et E e. —0.13 68 
293 3 5 8 102 4 1/2 0.456 et e e 1.48 

12229+2643 3716 3 5 6 76 2 1 0.280 e E e. 0.02 69 
151 3 5 6 TI 5 1/2 0.231 e E -- -0.07 

12234--0611 3566 3 5 8 66 2 1 0.200 0.071 0.038 —1.2 >2 70 
133 3 5 8 66 0 1/2 0.061 e E e 1.90 

J2234+0944 7535 3 5 5 72 2 1 0.197 ... E e. —0.17 71 
245 3 5 5 74 0 1/2 0.796 0.209 0.176 —0.1 >2 

12302+4442 10,211 3 5 7 108 3 1 0.764 et E -- —0.05 74 
236 3 5 7 108 0 1/2 0.710 et E e —0.03 

J23174-1439 13,942 3 5 6 303 3 2 0.345 e E et —0.09 73 
711 3 5 6 309 0 2/3 0.690 et E e 0.01 

J23224-2057 3088 3 5 0 59 1 2 0.255 et E e —0.25 74 
130 3 5 0 59 0 2/3 0.262 et E e. —0.13 

Notes. 


* The first line for each pulsar is from the narrowband analysis, and the second line is from the wideband analysis. 

Fit parameters: S — spin; A — astrometry; B — binary; DM — dispersion measure; FD — frequency dependence; J — jump (two numbers indicate wideband data 
with JUMPs/DMJUMP;). 
€ Weighted rms of epoch-averaged post-fit timing residuals, calculated using the procedure described in Appendix D of NG9. For sources with RN, the "Full" rms 
value includes the RN contribution, while the “White” rms does not. 
d RN parameters: Ад = amplitude of RN spectrum at = 1 уг! measured in pis yr 3, Vea = Spectral index; В = Bayes factor (“>2” indicates a Bayes factor larger 
than our threshold logio B > 2, but which could not be estimated using the Savage-Dickey ratio). See Equation (3) and Appendix С of NG9 for details. 
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Appendix C 


A complete set of timing residuals and DM variations for 
each pulsar in our data set can be found in Appendix C. 
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Figure 7. Narrowband and wideband timing residuals and DMX time series for J00234-0923. The wideband data set (and therefore the wideband residuals and DMX 
time series presented here) are only used for pipeline development and the mass determinations presented in Section 5.1.3. Top panel: narrowband and wideband 
DMX variations. Second panel: residual arrival times for all TOAs. Points are semitransparent; dark regions arise from the overlap of many points. Third panel: 
average residual arrival times. Bottom panel: all wideband TOA residuals. Receivers and back ends corresponding to each TOA are shown in the bottom panel’s 
legend. Dashed vertical line(s) denote the divide between back-end systems. 
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Figure 8. Narrowband and wideband timing residuals and DMX time series for 10030--0451. The wideband data set (and therefore the wideband residuals and DMX 
time series presented here) are only used for pipeline development and the mass determinations presented in Section 5.1.3. Top panel: narrowband and wideband 
DMX variations. Second panel: residual arrival times for all TOAs. Points are semitransparent; dark regions arise from the overlap of many points. Third panel: 
average residual arrival times. Fourth panel: whitened average narrowband TOA residuals. Fifth panel: all wideband TOA residuals. Bottom panel: all whitened 
wideband TOA residuals. Receivers and back ends corresponding to each TOA are shown in the bottom panel’s legend. Dashed vertical line(s) denote the divide 


between back-end systems. 
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Figure 9. Narrowband and wideband timing residuals апа DMX time series for J03404-4130. See Figure 7 for details. 
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Figure 10. Narrowband and wideband timing residuals and DMX time series for J0406+3039. See Figure 7 for details. 
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Figure 11. Narrowband and wideband timing residuals and DMX time series for J0437—4715. See Figure 8 for details. 
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Figure 12. Narrowband and wideband timing residuals and DMX time series for J05094-0856. See Figure 7 for details. 


Ј0557+1551 PUPPI 


NB & WB DMX 
[1073 рс ст“ 3] 
h о 


ы 
= 


All NB TOA 
Residual [us] 
о 


Avg. NB TOA 
Residual [us] 
b © 


| 
щл 
ө © 


All WB TOA 
Residual [us] 
© 


ч * PUPPIL * PUPPIS 


2004 2006 2008 2010 2012 2014 2016 2018 2020 
Year 


Figure 13. Narrowband and wideband timing residuals and DMX time series for J05574-1551. See Figure 7 for details. 


30 


THE ASTROPHYSICAL JOURNAL LETTERS, 951:L9 (78pp), 2023 July 1 Agazie et al. 


J0605+3757 GUPPI 


Рафи " 


All NB TOA NB & WB DMX 
Residual [us] [1073 pe ст“ 3] 


Avg. NB TOA 
Residual [us] 


All WB TOA 
Residual [us] 


+ СОРРИ. + GUPPI 820 MHz 


2004 2006 2008 2010 2012 2014 2016 2018 2020 
Year 


Figure 14. Narrowband and wideband timing residuals and DMX time series for Ј0605+3757. See Figure 7 for details. 
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Figure 15. Narrowband and wideband timing residuals and DMX time series for Ј0610–2100. See Figure 8 for details. 
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Figure 16. Narrowband and wideband timing residuals and DMX time series for Ј0613–0200. See Figure 8 for details. 
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Figure 17. Narrowband and wideband timing residuals and DMX time series for J0614—3329. See Figure 7 for details. 
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Figure 18. Narrowband and wideband timing residuals and DMX time series for Ј0636+5128. See Figure 7 for details. 
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Figure 19. Narrowband and wideband timing residuals and DMX time series for J06454-5158. See Figure 7 for details. 
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Figure 20. Narrowband and wideband timing residuals and DMX time series for JO709+0458. See Figure 7 for details. 
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Figure 21. Narrowband and wideband timing residuals and DMX time series for J07404-6620. See Figure 7 for details. 
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Figure 22. Narrowband and wideband timing residuals and DMX time series for J0931—1902. See Figure 7 for details. 
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Figure 23. Narrowband and wideband timing residuals and DMX time series for J1012+5307. See Figure 8 for details. 
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Figure 24. Narrowband and wideband timing residuals and DMX time series for J1012—4235. See Figure 7 for details. 
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Figure 25. Narrowband and wideband timing residuals and DMX time series for J1022+1001. See Figure 7 for details. 
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Figure 26. Narrowband and wideband timing residuals and DMX time series for J1024—0719. See Figure 7 for details. 
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Figure 27. Narrowband and wideband timing residuals and DMX time series for J1125+7819. See Figure 7 for details. 
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Figure 28. Narrowband and wideband timing residuals and DMX time series for J1312+0051. See Figure 7 for details. 
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Figure 29. Narrowband and wideband timing residuals and DMX time series for J14534-1902. See Figure 7 for details. 
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Figure 30. Narrowband and wideband timing residuals and DMX time series for 1455-3330. See Figure 7 for details. 


41 


THE ASTROPHYSICAL JOURNAL LETTERS, 951:L9 (78pp), 2023 July 1 Agazie et al. 


J1600—3053 GASP GUPPI/YUPPI 


4 NB DMX 


0 к-к. C аила каканак T - 


i кыы 
i x %х ” 
-2 Salen ұсы, relier 2 


50 


NB & WB DMX 
[1073 рс ст“ 3] 


E I: : — 


All NB TOA 
Residual [us] 
© 


—50 


Avg. NB TOA 
Residual [us] 
[=] 


| 
л 


Whitened Avg. 
NB TOA 
Residual [us] 

di b © N л 
© Ш © Un © 


All WB TOA 
Residual [us] 
© 


—4 

az 
== 2 
э 
о 
59 

“ 2 
кш 
< g =: 

2004 2006 2008 2010 2012 2014 2016 2018 2020 
Уеаг 


Figure 31. Narrowband and wideband timing residuals and DMX time series for J1600—3053. See Figure 8 for details. 
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Figure 32. Narrowband and wideband timing residuals and DMX time series for J1614—2230. See Figure 8 for details. 
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Figure 33. Narrowband and wideband timing residuals and DMX time series for J1630+3734. See Figure 7 for details. 
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Figure 34. Narrowband and wideband timing residuals and DMX time series for J1640+2224. See Figure 7 for details. 
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Figure 35. Narrowband and wideband timing residuals and DMX time series for J1643—1224. See Figure 8 for details. 
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Figure 36. Narrowband and wideband timing residuals and DMX time series for J1705—1903. See Figure 8 for details. 
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Figure 37. Narrowband and wideband timing residuals and DMX time series for J1713+0747. See Figure 8 for details. 
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Figure 38. Narrowband and wideband timing residuals and DMX time series for J1719—1438. See Figure 7 for details. 
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Figure 39. Narrowband and wideband timing residuals and DMX time series for Ј1730–2304. See Figure 7 for details. 
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Figure 40. Narrowband and wideband timing residuals and DMX time series for Ј1738+0333. See Figure 7 for details. 
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Figure 41. Narrowband and wideband timing residuals and DMX time series for J1741+1351. See Figure 7 for details. 
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Figure 42. Narrowband and wideband timing residuals and DMX time series for J1744—1134. See Figure 8 for details. 


50 


Авале et al. 


THE ASTROPHYSICAL JOURNAL LETTERS, 951:L9 (78pp), 2023 July 1 


Ј1745+1017 РОРРГ 


NB & WB DMX 
[1073 рс ст“ 3] 


АПХВ ТОА 
Residual [Us] 


Avg. NB TOA 
Residual [us] 


| 
> 


bh = 
245 
yO 
oes 
& 2 
525 
= = 


All WB TOA 
Residual [us] 
© 


25 
эч 
53 
2% 
E 2 
= < 
28 
2004 2006 2008 2010 2012 2014 2016 2018 2020 
Үсаг 


Figure 43. Narrowband and wideband timing residuals and DMX time series for J1745+1017. See Figure 8 for details. 
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Figure 44. Narrowband and wideband timing residuals and DMX time series for J1747—4036. See Figure 8 for details. 
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Figure 45. Narrowband and wideband timing residuals and DMX time series for Ј1751–2857. See Figure 7 for details. 
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Figure 46. Narrowband and wideband timing residuals and DMX time series for J1802—2124. See Figure 8 for details. 
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Figure 47. Narrowband and wideband timing residuals and DMX time series for Ј1811–2405. See Figure 7 for details. 
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Figure 48. Narrowband and wideband timing residuals and DMX time series for 1832-0836. See Figure 7 for details. 
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Figure 49. Narrowband and wideband timing residuals and DMX time series for J1843-1113. See Figure 7 for details. 
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Figure 50. Narrowband and wideband timing residuals and DMX time series for Ј1853+1303. See Figure 8 for details. 
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Figure 51. Narrowband and wideband timing residuals and DMX time series for B1855+09. See Figure 8 for details. 
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Figure 52. Narrowband and wideband timing residuals and DMX time series for J1903+0327. See Figure 8 for details. 
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Figure 53. Narrowband and wideband timing residuals and DMX time series for Ј1909–3744. See Figure 8 for details. 
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Figure 54. Narrowband and wideband timing residuals and DMX time series for J19104-1256. See Figure 7 for details. 
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Figure 55. Narrowband and wideband timing residuals and DMX time series for J19114-1347. See Figure 7 for details. 
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Figure 56. Narrowband and wideband timing residuals and DMX time series for Ј1918–0642. See Figure 7 for details. 


J192342515 ASP PUPPI 


NB & WB DMX 
[103 рс ст“ 3] 


АПХВ ТОА 
Residual [us] 
| 
8 e 


= 


Avg. NB TOA 
Residual [us] 


All WB TOA 
Residual [us] 


^ ASP 430MHz ` PUPPI430MHz | ASPL . PUPPIL 


2004 2006 2008 2010 2012 2014 2016 2018 2020 
Year 


Figure 57. Narrowband and wideband timing residuals and DMX time series for J1923+2515. See Figure 7 for details. 
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Figure 58. Narrowband and wideband timing residuals and DMX time series for B1937+21. See Figure 8 for details. 
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Figure 59. Narrowband and wideband timing residuals and DMX time series for J19444-0907. See Figure 7 for details. 
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Figure 60. Narrowband and wideband timing residuals and DMX time series for J1946+3417. See Figure 8 for details. 
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Figure 61. Narrowband and wideband timing residuals and DMX time series for B1953+29. See Figure 8 for details. 
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Figure 62. Narrowband and wideband timing residuals and DMX time series for J2010-1323. See Figure 7 for details. 
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Figure 63. Narrowband and wideband timing residuals and DMX time series for J2017+0603. See Figure 7 for details. 
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Figure 64. Narrowband and wideband timing residuals and DMX time series for J2033+1734. See Figure 7 for details. 
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Figure 65. Narrowband and wideband timing residuals апа DMX time series for J2043+1711. See Figure 7 for details. 
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Figure 66. Narrowband and wideband timing residuals and DMX time series for J2124—3358. See Figure 7 for details. 
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Figure 67. Narrowband and wideband timing residuals and DMX time series for ]2145-0750. See Figure 8 for details. 
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Figure 68. Narrowband and wideband timing residuals and DMX time series for J2214+3000. See Figure 7 for details. 
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Figure 69. Narrowband and wideband timing residuals and DMX time series for J2229+2643. See Figure 7 for details. 
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Figure 70. Narrowband and wideband timing residuals and DMX time series for J2234+0611. See Figure 8 for details. 
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Figure 71. Narrowband and wideband timing residuals and DMX time series for J2234+0944. See Figure 7 for details. 
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Figure 72. Narrowband and wideband timing residuals and DMX time series for J2302+4442. See Figure 7 for details. 
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Figure 73. Narrowband and wideband timing residuals and DMX time series for J2317+1439. See Figure 7 for details. 
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Figure 74. Narrowband and wideband timing residuals and DMX time series for J2322+2057. See Figure 7 for details. 
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Appendix D 
Kopeikin Parameter Implementations 


For PSR J1713--0747, which is one of the most precisely 
timed pulsars, we use parameters that ultimately relate to the 
orientation of the binary orbit on the sky. These are often called 
the “Kopeikin parameters” after two papers (Kopei- 
kin 1995, 1996) that laid out the relevant expressions. 

We begin with a recap of parameters used in pulsar timing 
codes: 


1. x—projected semimajor axis of the pulsar’s orbit; 

2. w—longitude of periastron; 

3. d—distance to the pulsar (derived from the parallax 
parameter); 

4. џа—рторег motion in R.A., а cos ё; 

5. ив-ргорег motion in decl., 6; 

6. i—orbital inclination angle (generally has a quadrant 
ambiguity, as often only sini is accessible through 
timing; has different definitions); 

7. Q—longitude of the ascending node (angle between the 
orbital plane and a reference line on the sky; this has 
different definitions as discussed below). 


Kopeikin (1995) defines the AOP term as follows (Equation 
(17)): 


Алм = СТ 5ш — Aj, cosQ)R(u)coti 


— (A, cosQ + Aj, sin О) Q(u)csci], (81) 


where О(и) is an orbital “cosine” term often designated as C in 
timing codes, К(и) is an orbital “sine” term often designated as 
S in timing codes, и is the eccentric anomaly, and А, and А, 
are dot products of Earth’s position with unit vectors pointing 
east and north, respectively, on the sky at the position of the 
pulsar system. Kopeikin (1995) points out that this term can be 
broken into effects on x and w: 


coti 
d 


ер CSCI 2 
(988 = jinrinsic _ (Ap cosQ + A, sinQ). 
а lo Jo 


obs — intrinsic | 1+ (Aj, sinQ — А „соз О) | (В2) 


(B3) 


Kopeikin (1996; see also Arzoumanian et al. 1996) presents 
the excess Romer delay due to the effects of changing projected 
orbit due to proper motion: 


EM бту) = x(t — to)csci (u, сов О + ps sin О) 
с 
x (1 — ecosu)cos(w + А.(и)) 
+ x(t — to)coti(—p, sin Q + us cos) 
x (1 — ecosu)sin(w + A,(u)), 
(B4) 


where r, is the changing radius vector of the pulsar’s orbit 
about the binary center of mass and А.(и) is the corresponding 
true anomaly. The resulting equations for the changes in 
observables are 


ôi = (— ua Sin Q + џи; cos 0) (7 — to) (85) 
е = coti(—p, sin Q + ps cosQ)(t — to) (B6) 
x 
бш = све (u, cos О + uis sinQ)(t — to), (B7) 
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where t — tọ is the time elapsed since the reference epoch. 

Between one and four of these effects can be observed in 
many systems, and all can be parameterized using i and О. 
Typically X? values are computed for a grid of i and О values; 
only strong “Kopeikin effects” will break the natural fourfold 
degeneracy in the parameters and narrow down the possible 
values of i and О. 

Unfortunately, there are two different sign conventions in 
use to describe i and О, with literature and software packages 
(Tempo, Tempo2, and PINT) using both. Use of the incorrect 
sign convention when updating pulsar parameters will result in 
poor fits and possible settling in a local rather than global 
minimum in the і-О space. 

To forestall such errors, we summarize the two conventions 
in use: 


1. IAU convention: 

(a) i— 0? means that the orbital angular momentum 
vector points toward Earth, and i= 180? means that 
the orbital angular momentum vector points away 
from Earth. 

(b) О is 0° toward the north and increases counter- 
clockwise on the sky; it is measured “north 
through east." 

2. Damour & Taylor (1992) (DT92) convention: 

(a) i= 180? means that the orbital angular momentum 
vector points toward Earth, and i — 0? means that the 
orbital angular momentum vector points away from 
Earth. 

(b) О is 0? toward the east and increases clockwise on the 
sky; it is measured "east through north." 


These conventions are related by the following: 


1. ipro? = 180° — йди; 
sin ipro2 = SiN iray; COS утог = — COS А 
o . 
2. Ортог = 90° — Ордо; 
sin Ортог = cosQyay; cos ртог = sin Орду. 


The ТАЏ convention has been used in Splaver et al. (2005) 
and NANOGrav Tempo-derived papers such as Fonseca et al. 
(2016) and the 12.5yr data release (Alam et al. 2021a). 
The DT92 convention was used to derive the equations in 
Kopeikin (1995, 1996) and has been used for measurements in 
van Straten et al. (2001) and Stairs et al. (2004). We include a 
diagram illustrating both conventions in Figure 6, which is 
adapted from Splaver et al. (2005). 

The conventions used in the various pulsar timing codes are 
as follows: 


1. Tempo DDK model (van Straten 2003): The code used to 
read and write parfiles uses the IAU convention. 
Internally to the code, the input i (“КІМ”) and О 
(“КОМ”) values are immediately transformed to 
the DT92 convention, and the Kopeikin equations are 
used directly in that convention. For each TOA, this code 
first computes the locally adjusted x and w based on 
Equations (B6) and (B7) followed by the AOP 
corrections to x and w based on Equations (B2) and 
(B3). It then proceeds through the standard computation 
of orbital delays, including Shapiro delay based on i, 
calculation of all relevant derivatives, and parameter 
adjustment. The fit results are transformed back to the 
IAU convention before output. Note that Tempo requires 
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а раге line 77К96 1”” to ensure that the proper- 
motion corrections are done. 

2. Tempo2 DDK model: The code used to read and write 

parfiles uses the DT92 convention for KIN and KOM. 

The delay /derivative code closely follows the logic of the 

Tempo DDK code, including the use of the DT92 
convention in the Kopeikin equations. We have recently 
debugged this routine and confirmed that it returns 
equivalent results to the tempo DDK model. 

3. Tempo2 T2 model: The code used to read and write 
parfiles uses the IAU convention for KIN and KOM. This 
model provides a superset of multiple binary models and 
uses Equation (B1) directly rather than breaking it down. 
The fitting code uses an IAU convention implementation 
of the Kopeikin equations. We note that only the base 
values of i, x, and w, not values adjusted via 
Equations (B5)-(B7), are used in determining the AOP 
terms and also the Romer delay. This introduces small 
discrepancies between its results and those of DDK, with 
the T2 model parameters being incorrect by up to roughly 
0.20, according to our simulations. 

4. PINT DDK model: The code used to read and write 
parfiles uses the DT92 convention for KIN and KOM. 
Internally the code uses the equations from Kopeikin 
(1995, 1996) directly. The PINT code does apply 
corrections based on Equations (B5)-(B7) before com- 
puting the AOP terms. 


In this paper, we present astrometric results in ecliptic 
coordinates rather than equatorial Tempo and PINT are 
written to fit for Q relative to ecliptic north (Tempo) or east 
(PINT), rather than equatorial north/east in this case, and use 
the proper motions in ecliptic coordinates (ил, ив) to calculate 
the relevant parameters. 
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